In [1]:
# Parameters
MCC = 1.5
BVALL = 0.84
SWARM_PATH = "/home/kamran.karimi1/Project/git/Swarm/dataset/YuhaDesert/EMC.csv"
In [2]:
#--- extract the catalog
import pandas as pd
#import matplotlib; matplotlib.use('agg')
import sys
from matplotlib import rc as rcc
import matplotlib
import matplotlib.pylab as plt
import numpy as np
from matplotlib import patches
import warnings
import matplotlib.ticker
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import datetime
#import path
from math import *
import matplotlib.dates as mdates
from matplotlib import dates
import matplotlib.dates as md
import matplotlib.cm as cm
from matplotlib.font_manager import FontProperties
import itertools
import scipy
from scipy.ndimage import gaussian_filter
import math
import pdb
def DrawFrame(ax, (alpha_xl,alpha_xr),(alpha_yb,alpha_yt),linewidth,LOG_X=None,LOG_Y=None):
[xlo,xhi,ylo,yhi] = ax.axis()
if LOG_X:
[xlo,xhi,junk,junk] = np.log10(ax.axis())
if LOG_Y:
[junk,junk,ylo,yhi] = np.log10(ax.axis())
lx = xhi - xlo
ly = yhi - ylo
xy = [xlo - alpha_xl * lx, ylo - alpha_yb * ly]
height = ly*(1+alpha_yb+alpha_yt)
width = lx*(1+alpha_xl+alpha_xr)
xy_end=[xy[0]+width,xy[1]+height]
if LOG_X:
xy[0] = 10 ** xy[0]
xy_end[0] = 10 ** xy_end[0]
if LOG_Y:
xy[1] = 10 ** xy[1]
xy_end[1] = 10 ** xy_end[1]
ax.add_patch( patches.Rectangle(xy=xy, width=xy_end[0]-xy[0],
height=xy_end[1]-xy[1], linewidth=linewidth,
clip_on=False,facecolor=None,edgecolor='black',fill=None) )
#--- add a new time attribute
def ConvertTime( df_in ):
df=df_in.copy()
df.insert(0,'date',pd.to_datetime(swarm[['year', 'month', 'day', 'hour', 'minute', 'second']]))
df.drop(['year', 'month', 'day', 'hour', 'minute', 'second'],axis=1,inplace=True)
return df
#--- set path
#SWARM_PATH = './dataset/YuhaDesert/EMC.csv' #--- comment if passed by arg
DIR_OUTPT = '.' #--- uncomment if passed by arg
#DIR_OUTPT = 'swarmYuhaDesertFirstPeriod0/Run0'
DIR_OUTPT_figs = '.'
#--- store
swarm = pd.read_csv( SWARM_PATH, sep = ' ' ) #--- parse data
swarm = ConvertTime( swarm ) #--- add new column 'date'
#--- sort based on time
swarm.sort_values(by=['date'],inplace=True)
#--- reindex
swarm.reset_index(inplace=True,drop=True)
#--- rename cols
swarm = swarm.rename(index=str,columns={'lat':'latitude','lon':'longitude','mag':'magnitude'})
swarm.sort_values(by='magnitude',ascending=None).head()
#pdb.set_trace()
Out[2]:
In [3]:
mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings
rcc('text', usetex=True)
font = {'size' : 20}
matplotlib.rc('font', **font)
In [4]:
#--- remove nan
#swarm = swarm.dropna(subset=['magnitude'])
swarm.info()
In [5]:
#--- plot spatial map
#--- fault map california: temporal evolution of events
def DataFrameSubSet( df, column, (xmin,xmax) ):
return df[ ( df[column] >= xmin ) &
( df[column] < xmax ) ]
#--- subset of data
swarm_lohi = swarm.copy()#DataFrameSubSet( swarm,
# 'date',
# ( pd.to_datetime('2010-04-04'), pd.to_datetime('2010-06-26') ) )
swarm_lohi.plot.scatter('longitude','latitude',
s=3**(swarm_lohi['magnitude']),
c='date',cmap='jet',alpha=0.1) #--- plot
#plt.savefig('%s/map.png'%DIR_OUTPT_figs,dpi=150,bbox_inches='tight')
#plt.figure(figsize=(6,6)).gca(projection='3d')
#plt.xlabel('Long')
#plt.ylabel('Lat')
#plt.scatter(swarm_lohi['longitude'],
# swarm_lohi['latitude'],
# swarm_lohi['depth']) #s=10*swarm['magnitude'],cmap='jet') #--- plot
plt.show()
plot spatial map¶
In [6]:
from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
import matplotlib.font_manager as fm
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,
myTlohi.hi ) )
exponent = 2.5
dff = swarm_lohi.copy()
# (injection_period, xlo, xhi ) = (0, dff['date'].min(),
# pd.to_datetime('2019-07-6 00:00:00'))
(injection_period, xlo, xhi ) = (1, dff['date'].min(),
dff['date'].max()
)
#--- plot
fig=plt.figure(figsize=(5,5))
ax=fig.add_subplot(111)
# ax.set_xlabel('Longitude')
# ax.set_ylabel('Latitude')
#ax.tick_params(labelsize='small')
#--- xmax, xmin
Xlo = swarm_lohi['longitude'].min()
Xhi = swarm_lohi['longitude'].max()
Ylo = swarm_lohi['latitude'].min()
Yhi = swarm_lohi['latitude'].max()
dx=0.5*max(Xhi-Xlo,Yhi-Ylo)*.75
xc=0.5*(Xlo+Xhi)
yc=0.5*(Ylo+Yhi)
ax.set_xlim(xc-dx,xc+dx)
ax.set_ylim(yc-dx,yc+dx)
#--- set tick labels
xx={0:xc,1:yc}
xlist=np.linspace(xx[0]-dx,xx[0]+dx,6,endpoint=True)
ax.set_xticks(xlist)
xlist=[r'$%4.2f^{\circ}$'%i for i in xlist]
#
ylist=np.linspace(xx[1]-dx,xx[1]+dx,6,endpoint=True)
ax.set_yticks(ylist)
ylist=[r'$%4.2f^{\circ}$'%i for i in ylist]
#ax.tick_params(labelsize=12)
ax.tick_params(axis='y',labelright=False,labelleft=True,right=True,left=True)
ax.tick_params(axis='x',top=True)#,left=True)
ax.set_xticklabels(xlist,{'fontsize':14})
ax.set_yticklabels(ylist,{'fontsize':14})
#--- filter
dff = dff[(dff['date'] > xlo) &\
(dff['date'] < xhi)] #--- inj
#--- plot
# ax.set_title(r'$%s$'df_lohi.iloc[itime]['date'])#,labelsize='small')
line=ax.scatter(dff['longitude'],dff['latitude'],
s=exponent**(dff['magnitude']),
# label=df_lohi.iloc[itime]['date'],
alpha=0.1,animated=True,color='C0') #--- plot
#ax.legend()
#--- plt grid
plt.grid()
# #--- mmax
# m_max = dff['magnitude'].max()
# x0=dff['longitude'][dff['magnitude']==m_max].iloc[0]
# y0=dff['latitude'][dff['magnitude']==m_max].iloc[0]
# plt.scatter(x0,y0,
# s=exponent**m_max,
# marker='*',zorder=3,
# facecolors='yellow',color='black',
# alpha=1.0)
# # #--- 2nd max
# m_max=dff.sort_values(by='magnitude',ascending=False)['magnitude'].iloc[1]
# x0=dff['longitude'][dff['magnitude']==m_max].iloc[0]
# y0=dff['latitude'][dff['magnitude']==m_max].iloc[0]
# plt.scatter(x0,y0,
# s=exponent**m_max,
# marker='*',zorder=3,
# facecolors='yellow',color='black',
# alpha=1.0)
fontprops = fm.FontProperties(size=14)
xx=1/111.0
yy = 5 #---km
scalebar = AnchoredSizeBar(ax.transData,
yy*xx, r'$%s~km$'%yy, 'upper right',label_top=True,
# pad=0.1,
fill_bar='white',
frameon=False,#,
size_vertical=yy*xx/10,
fontproperties=fontprops)
ax.add_artist(scalebar)
DrawFrame(ax, (0.17,0.09),(0.15,0.06),0.01)
fig.savefig('%s/spatialMapYuha.png'%(DIR_OUTPT_figs),dpi=150,bbox_inches='tight')
plt.show()
In [7]:
#--- plot timeseries
class Tlohi:
def __init__(self,lo,hi):
self.lo = lo
self.hi = hi
myTlohi = Tlohi( pd.to_datetime('2010-04-04 00:00:00'),
pd.to_datetime('2010-04-14'))
swarm_lohi = DataFrameSubSet(swarm,'date',(myTlohi.lo,myTlohi.hi))
fig=plt.figure(figsize=(8,4))
ax=fig.add_subplot(111)
plt.xlabel('Time')
plt.ylabel('M')
plt.ylim(-2,swarm['magnitude'].max())
plt.xlim(myTlohi.lo, myTlohi.hi) #pd.to_datetime('2010-04-04'),pd.to_datetime('2010-06-26'))
#plt.xlim(swarm['date'].min(),pd.to_datetime('2010-04-24'))
#plt.xlim(myTlohi.lo,myTlohi.hi)
plt.scatter(swarm_lohi['date'],swarm_lohi['magnitude'],
s=np.exp(swarm_lohi['magnitude']),
alpha=.1) #0.02)
# plt.plot([myTlohi.lo,myTlohi.lo],
# [-2,7],'r-')
# plt.plot([myTlohi.hi,myTlohi.hi],
# [-2,7],'r-')
fig.autofmt_xdate()
plt.savefig('%s/timeSeries.png'%DIR_OUTPT_figs,dpi=150,bbox_inches='tight')
plt.show()
In [8]:
#--- freq-magnitude relation
import sys
def histogramACCUMLTD( slist ):
assert type( slist ) == type( [] ), 'arg must be a list. a %s is given!' %( type( slist ) )
d = {}
for item in slist:
try:
d[ item ] += 1
except:
d[ item ] = 1
keys = d.keys()
keys.sort()
cdf = 0.0
xi = min( slist ) - 1.0e-6
xf = max( slist ) + 1.0e-6
npoin = len( slist )
adict = {}
for ikey, index in zip( keys, xrange( sys.maxint ) ):
adict[ index ] = [ xi, ikey, cdf ]
cdf += 1.0 * d[ ikey ] # / npoin
xi = ikey
adict[ index + 1 ] = [ xi, xf, cdf ]
return adict
#--- set min/max time to avoid temporal incompletenesss issue
#myTlohi.hi = pd.to_datetime('2010-04-08 00:00:00')
swarm_copy = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,
myTlohi.hi ) )
#--- accumulated histogram
N = len(swarm_copy['magnitude'])
slist=np.array(swarm_copy['magnitude'])
slist.sort()
d = histogramACCUMLTD( slist.tolist() )
keys=d.keys()
keys.sort()
#--- plot
fig= plt.figure(figsize=(3,3))#,dpi=150)
ax = fig.add_subplot(111)
ax.axis([-2,6,1,1e6])
ax.set_yscale('log')
#--- add major xticks
xmin=np.ceil(ax.axis()[0])
xmax=np.floor(ax.axis()[1])
nbin = xmax - xmin
ax.set_xticks(np.linspace(ax.axis()[0],ax.axis()[1],int(nbin)+1))
#--- add major yticks
ymin=np.ceil(np.log10(ax.axis()[2]))
ymax=np.floor(np.log10(ax.axis()[3]))
nbin = ymax - ymin
ax.set_yticks(np.logspace(ymin,ymax,int(nbin)+1))
#--- put minor bins
#locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
#ax.xaxis.set_minor_locator(locmin)
#ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
#--- put minor bins y
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
ax.tick_params(axis='y',left=True, right=True,which='both',labelleft=False, labelright=True)
ax.tick_params(axis='x',bottom=True, top=True,which='both')
xmin=swarm_copy['magnitude'].min()
xmax=swarm_copy['magnitude'].max()
dx=0.1
junk = ax.hist( swarm_copy['magnitude'],
bins=int((xmax-xmin)/dx),
label='histogram',color='silver') #--- histogram
xx=[];yy=[]
for ikey in keys:
xx.append(d[ikey][0])
yy.append(N-d[ikey][2])
ax.plot(xx,yy,
linestyle='-', drawstyle='steps-post',color='black',
linewidth=1.0) #--- accumulated
b=bval# 0.84#0.95
ax.plot([0,6],[1e6, 1e6* 10**(-b*(6-0))],'r-.',linewidth=1)
DrawFrame(ax, (0.08,0.2),(0.32,0.1),0.01,LOG_Y=True)
ax.set_title(r'%s~%s-%s~%s, b=%s'%(myTlohi.lo.month_name(),myTlohi.lo.day,
myTlohi.hi.month_name(),myTlohi.hi.day,b),fontsize='xx-small')
plt.savefig('%s/gr.png'%DIR_OUTPT_figs,dpi=150,bbox_inches='tight')#,transparent=True)
plt.show()
In [9]:
#--- estimate b by varying the mc value
from math import *
import scipy.stats
#--- set min/max time to avoid temporal incompletenesss issue
swarm_copy = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,
myTlohi.hi ) )
#--- plot distributions (Nm)
fig = plt.figure(figsize=(12,4))
ax=fig.add_subplot(131)
ax.set_xlabel('$m$')
ax.set_ylabel('$N(mag.>m)$')
xmin = -3 #np.floor(swarm_copy['magnitude'].min())
xmax = 5 #np.ceil(swarm_copy['magnitude'].min())
ax.set_xticks(np.linspace(xmin,xmax,int(xmax-xmin)/1+1))
ax.set_xlim(xmin,xmax)
#plt.ylim(1,1e5)
ax.set_yscale('log')
for ikey in keys:
ax.plot([d[ikey][0],d[ikey][1]],
[N-d[ikey][2],N-d[ikey][2]],
'-',color='black') #--- accumulated
junk = ax.hist( swarm_copy['magnitude'],
bins=64,
label='histogram') #--- histogram
#--- list of magnitudes
mm_list = swarm_copy['magnitude'].to_list()
mm_list = np.sort(mm_list)
#--- b vs mc
ax2=fig.add_subplot(132)
ax2.set_xticks(np.linspace(-1.0,3.0,5))
#ax4 = ax2.twinx()
#ax4.set_yscale('log')
#ax4.set_ylim(1e-16,)
#plt.xlim(0.0,1.5)
ax2.set_xlabel('$m_c$')
ax2.set_ylabel('$b$')
ax2.set_ylim(0.4,1.2)
#ax4.set_ylabel(r'$p_{val}$',labelpad=-70)
#--- bvalue bhatacharya eq. 5
DM = 0.1
n = 51
#---
MC = []
B = []
Err_B = []
pval = []
for mc_new in [-1.05+i*0.1 for i in xrange(n)]: #--- m values are quantized (+- 0.1)!
mlist = swarm_copy[ swarm_copy['magnitude']-mc_new>=0.0]['magnitude']
m_mean=mlist.mean()
assert m_mean-mc_new > 0.0, 'm_mean=%s,mc=%s'%( m_mean, mc )
k_mean = (m_mean-mc_new)/DM
b_value = log(np.exp(1),10)*log(1+DM/(m_mean-mc_new))/DM
nsamp = len(mlist)
# err_b = 1.96*(b_value/nsamp/log(10,exp(1.0)))**0.5
err_b = 2*(k_mean*(k_mean+1)/nsamp)**0.5/log(10,exp(1.0)) #--- Guttorp and Hopkins (1986) eq.(4)
#--- bvalue bhatacharya eq. 6
#--- accumulated histogram
# slist = mm_list[ mm_list>=mc_new ]
# slist.sort()
# NN = len(slist)
# cdf = histogramACCUMLTD( slist.tolist() )
# # pdb.set_trace()
# #--- bvalue bhatacharya eq. 6
# n_cdf=len(cdf)-1
# bi = []
# for i in xrange(n_cdf):
# mi = cdf[ i ][ 0 ]
# Ni = NN - cdf[ i ][ 2 ]
# bij = []
# for j in xrange(n_cdf):
# mj = cdf[ j ][ 0 ]
# Nj = NN - cdf[ j ][ 2 ]
# if i == j: continue
# # print i,j,mi,mj
# # try:
# slope = ( log( Nj, 10 ) - log( Ni, 10 ) ) / ( mj - mi )
# assert not np.isnan(slope)
# bij.append( slope )
# # except:
# # continue
# bi.append( np.median( bij ) )
# bval_2nd = -np.median( bi )
# err_b2nd = np.std(bi)/(len(bi))**0.5
#--- plot b vs mc
ax2.errorbar([mc_new],[b_value],yerr=err_b,marker='o',color='black')
# ax2.errorbar([mc_new],[bval_2nd],yerr=err_b2nd, marker='s',color='red')
#---
MC.append(mc_new)
B.append(b_value)
Err_B.append(err_b)
#--- p-value (by hand)
# (dmax,pval2nd) = scipy.stats.kstest( mlist, lambda x: 1.0-10**(-b_value*(x-mc_new)))
# ax4.plot([mc_new],[pval2nd],marker='o',color='red')
# pval.append(pval2nd)
#--- choose the range of mc
mc_range=[3,3]
index_tf = [i and j for i, j in zip(mc_range[0]<=np.array(MC),np.array(MC)<mc_range[1])]
bvall = np.array(B)[index_tf]
Err_B = np.array(Err_B)[index_tf]
#--- b vs mc: vertical lines
#ax2.plot([mc_range[0],mc_range[0]],[np.min(B),np.max(B)],'--r')
#ax2.plot([mc_range[1],mc_range[1]],[np.min(B),np.max(B)],'--r')
#--- N vs m: vertical line
#ax.plot([mc_range[0]+DM/2,mc_range[0]+DM/2],
# [1,len(swarm_copy['magnitude'])],'--r')
#--- n*10*bm (rescaled)
bval=1.3 #0.84#0.95 #np.mean(bvall)
#bval2=0.9
ax3=fig.add_subplot(133)
ax3.set_xlabel('m')
ax3.set_ylabel(r'$N\times10^{%s\times m}$'%bval)
#ax3.set_ylabel(r'$N\times10^{bm}$')
ax3.yaxis.set_label_position("right")
ax3.yaxis.tick_right()
ax3.set_yscale('log')
ax3.set_xticks(np.linspace(xmin,xmax,int(xmax-xmin)/1+1))
ax3.set_xlim(xmin,xmax)
ax3.set_ylim(1e3,1e6)
#ax3.set_xticks(np.linspace(ax3.axis()[0],ax3.axis()[1],4))
for ikey in keys: #--- rescaled
c = 10**(bval*d[ikey][0])
ax3.plot([d[ikey][0],d[ikey][1]],
[c*(N-d[ikey][2])/100,c*(N-d[ikey][2])/100],'-',color='black',
) #--- accumulated
# c = 10**(bval2*d[ikey][0])
# ax3.plot([d[ikey][0],d[ikey][1]],
# [c*(N-d[ikey][2]),c*(N-d[ikey][2])],'-',color='red',
# ) #--- accumulated
ax3.plot([mc_range[0]+DM/2,mc_range[0]+DM/2], #--- vertical line
[len(swarm_copy[swarm_copy['magnitude']>mc_range[0]])*10**(bval*mc_range[0])/10,
len(swarm_copy[swarm_copy['magnitude']>mc_range[0]])*10**(bval*mc_range[0])],'--r')
#plt.plot(MC, pval,'-o',color='black')
fig.tight_layout(w_pad=-0.0)
#fig.savefig('%s/b.png'%DIR_OUTPT_figs,dpi=75)
plt.show()
In [10]:
# #--- estimate b and mc
# #--- set min/max time to avoid temporal incompletenesss issue
# swarm_copy = swarm.copy()#DataFrameSubSet( swarm,
# # 'date',
# # ( myTlohi.lo,
# # myTlohi.hi ) )
# #--- power law fit clauset
# n_samp = 1000
# n_chunk = len( swarm_copy ) / n_samp
# i = 0
# fig = plt.figure(figsize=(8,4))
# ax = fig.add_subplot(111)
# ax.set_ylabel('m')
# ax2=ax.twinx()
# ax.scatter(swarm_copy['date'],swarm_copy['magnitude'],
# alpha=0.1)
# fig2 = plt.figure(figsize=(8,4))
# ax21 = fig2.add_subplot(111)
# ax21.set_ylabel('m')
# ax21.scatter(swarm_copy['date'],swarm_copy['magnitude'],
# alpha=0.1)
# MC = []
# B = []
# for j in xrange(n_chunk):
# Mlist = 10**(swarm_copy.iloc[i:i+n_samp]['magnitude'])
# fit = plfit.plfit( np.array(Mlist) ) #--- clauset
# mc = log(fit[1],10) #--- completeness mag.
# bval = fit[0]-1
# MC.append(mc)
# B.append(bval)
# # print mc, bval
# ax2.plot([swarm_copy.iloc[i:i+n_samp]['date'].min(),swarm_copy.iloc[i:i+n_samp]['date'].max()],
# [bval,bval],'-',color='red')
# ax21.plot([swarm_copy.iloc[i:i+n_samp]['date'].min(),swarm_copy.iloc[i:i+n_samp]['date'].max()],
# [mc,mc],'-',color='red')
# i += n_samp
# Mlist = 10**(swarm_copy.iloc[i:i+len(swarm_copy)%n_samp]['magnitude'])
# fit = plfit.plfit( np.array(Mlist) ) #--- clauset
# mc = log(fit[1],10) #--- completeness mag.
# bval = fit[0]-1
# MC.append(mc)
# B.append(bval)
# #---- plot bval
# ax2.plot([swarm_copy.iloc[i:i+n_samp]['date'].min(),swarm_copy.iloc[i:i+n_samp]['date'].max()],
# [bval,bval],'-',color='red')
# ax2.tick_params(colors='red')
# ax2.set_ylabel('b',color='red')
# ax.set_xlabel('Date')
# #print mc, bval
# #fig.savefig('%s/b.png'%DIR_OUTPT_figs,dpi=75)
# #--- plot mc
# ax21.plot([swarm_copy.iloc[i:i+n_samp]['date'].min(),swarm_copy.iloc[i:i+n_samp]['date'].max()],
# [mc,mc],'-',color='red',label='mc')
# plt.legend()
# ax21.set_xlabel('Date')
# #n, bins, patches = plt.hist(B, 50)
# #--- error estimate
# #[ error_b, error_mc, ntail ] = [ 0.0, 0.0, 0.0 ]
# #[error_b, error_mc, ntail] = plvar.plvar(np.array(Mlist), 'silent')
# #error_mc = log(1.0+error_mc/fit[1],10) #--- log(1+dM/M)
# #plt.savefig('%s/mc.png'%DIR_OUTPT_figs,dpi=75)
In [11]:
# #--- goodness of the fit
# n_samp = 1000
# n_chunk = len( swarm_copy ) / n_samp
# ncols=4
# nrows = n_chunk / ncols + int(np.ceil(1.0*( len(swarm_copy) % n_samp) / n_samp))
# i = 0
# plt.figure(figsize=(ncols*4,nrows*4))
# for j in xrange(n_chunk):
# Mlist = 10**(swarm_copy.iloc[i:i+n_samp]['magnitude'])
# fit = plfit.plfit( np.array(Mlist) ) #--- clauset
# mc = log(fit[1],10) #--- completeness mag.
# bval = fit[0]-1
# # [error_b, error_mc, ntail] = plvar.plvar(Mlist, 'silent') #--- error
# plt.subplot(nrows,ncols,j+1)
# Mlist=np.log10(Mlist)
# nc = len(Mlist[Mlist>=mc])
# error_b = (bval / nc / log(10.0) )**0.5
# hist,bins = np.histogram(Mlist,bins=1024)
# plt.plot(bins[:-1],-hist.cumsum()+len(Mlist))#,marker='o',markersize=6,linestyle='None')
# xc=0.5*(bins[:-1]+bins[1:])
# plt.plot(xc, len(Mlist[Mlist>=mc])*10**(-bval*(xc-mc)),'r',label='$b=%2.1f\pm%3.2f$'%(bval,error_b))
# plt.plot([mc,mc],[1,len(Mlist[Mlist>=mc])],label='$mc=%2.1f$'%(mc)) # len(Mlist[Mlist>=mc])*10**(-bval*(xc-mc)),'r')
# plt.yscale('log')
# plt.xlim(swarm_copy['magnitude'].min(),swarm_copy['magnitude'].max())
# plt.ylim(1,n_samp*10)
# if j == 0:
# plt.ylabel('N(size>m)')
# plt.xlabel('m')
# plt.legend()
# i += n_samp
# Mlist = 10**(swarm_copy.iloc[i:i+len(swarm_copy)%n_samp]['magnitude'])
# fit = plfit.plfit( np.array(Mlist) ) #--- clauset
# mc = log(fit[1],10) #--- completeness mag.
# bval = fit[0]-1
# Mlist=np.log10(Mlist)
# nc = len(Mlist[Mlist>=mc])
# error_b = (bval / nc / log(10.0) )**0.5
# hist,bins = np.histogram(Mlist,bins=1024)
# plt.subplot(nrows,ncols,j+2)
# plt.plot(bins[:-1],-hist.cumsum()+len(Mlist))#,marker='o',markersize=6,linestyle='None')
# xc=0.5*(bins[:-1]+bins[1:])
# plt.plot(xc, len(Mlist[Mlist>=mc])*10**(-bval*(xc-mc)),'r',label='$b=%2.1f\pm%3.2f$'%(bval,error_b))
# plt.plot([mc,mc],[1,len(Mlist[Mlist>=mc])],label='$m_c=%2.1f$'%(mc)) # len(Mlist[Mlist>=mc])*10**(-bval*(xc-mc)),'r')
# plt.yscale('log')
# plt.xlim(swarm_copy['magnitude'].min(),swarm_copy['magnitude'].max())
# plt.ylim(1,n_samp)
# plt.legend()
# #plt.savefig('%s/b_mult.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight')
In [12]:
#--- plot complete catalog
def ConvertDailyRate(hist, bin_edges ):
#---convert data to daily rate
t0 = pd.to_datetime( bin_edges[ 0 ] )
t1 = pd.to_datetime( bin_edges[ 1 ] )
delta_t = ( t1 - t0 ).total_seconds() / ( 60 * 60 * 24 )
hist *= ( bin_edges[ 1 ] - bin_edges[ 0 ] ) / delta_t
def ActivityRate( swarm ):
nbins = int( (swarm['date'].max()-swarm['date'].min()).days ) #--- number of bins
tmax = swarm['date'].max().value #--- min/max
tmin = swarm['date'].min().value
hist, bin_edges = np.histogram(swarm['date'].apply(lambda x:x.value),
bins=np.linspace(tmin,tmax,nbins+1,endpoint=True),density=True) #--- probability dist.
hist *= len( swarm['date'] ) #--- int(hist).dt=n
cumm_number = np.cumsum(hist)*(bin_edges[1]-bin_edges[0]) #--- accumulated number
ConvertDailyRate( hist, bin_edges ) #--- daily activity
return bin_edges, hist, cumm_number
#---------------------------------------------------------------------------------
#-----------------
#-----------------
#-----------------
#---------------------------------------------------------------------------------
#--- completeness
mc = 1.0
#--- t0<t<t1
#--- set min/max time to avoid temporal incompletenesss issue
swarm_tmp = DataFrameSubSet( swarm,
'date',
( swarm['date'].min(),
pd.to_datetime('2010-07-02') ) )
#--- m > mc
swarm_lohi = swarm_tmp.copy() #DataFrameSubSet( swarm_tmp,
# 'magnitude',
#( mc, sys.maxint ) )
#--- spatial map
#swarm_lohi.plot.scatter('longitude','latitude',
# s=3**(swarm_lohi['magnitude']),
# c='date',cmap='jet',
# alpha = 0.4) #--- plot
#--- temporal map
fig = plt.figure(figsize=(8,3))#,dpi=150)
ax1=fig.add_subplot(111)
#plt.xlabel('Time')
#plt.ylabel('M')
ax1.set_ylim(-1,6)
ax1.set_xlim(swarm_tmp['date'].min(),swarm_tmp['date'].max())
ax1.tick_params(axis='both',labelsize=18)
ax1.tick_params(axis='x',rotation=30)
ax1.scatter(swarm_lohi['date'],swarm_lohi['magnitude'],
s=2*np.exp(swarm_lohi['magnitude']),
alpha=0.04,color='black')
DrawFrame(ax1, (0.11,0.07),(0.32,0.1),0.01)
#--- activity rate
bin_edges, hist, cumm_number = ActivityRate( swarm_lohi )
#--- plot
ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis
ax2.plot(pd.to_datetime(bin_edges[:-1]),hist,'r-',linewidth=1)
ax2.tick_params(axis='y',labelsize=18,colors='red')
ax2.set_yscale('log')
ax2.set_ylim(1e2,1e4)
ax2.xaxis.set_major_formatter(dates.DateFormatter('%b %y'))
#--- add major xticks
xmin=ax2.axis()[0]
xmax=ax2.axis()[1]
nbin = 3
ax2.set_xticks(np.linspace(xmin,xmax,int(nbin)+1))
#fig.autofmt_xdate()
#fig.savefig('%s/timeSeriesYuhaDesert.png'%DIR_OUTPT_figs,dpi=150,bbox_inches='tight')
plt.show()
In [13]:
#--- evaluate fractal dimension
import geopy.distance
from math import *
import random as rnd
class bins:
def __init__(self, nBin, xlo, xhi, ydim = 1, err = None):
self.lo = xlo - 1.0e-10
self.hi = xhi + 1.0e-10
self.dx = (self.hi-self.lo)/nBin
self.xlist = [0.0 for i in xrange(nBin)]
self.kounter = [0 for i in xrange(nBin)]
self.ylist = [[0.0 for j in xrange(ydim)] for i in xrange(nBin)]
self.nBin = nBin
self.err = err
self.max_y = [[-sys.maxint for j in xrange(ydim)] for i in xrange(nBin) ]
self.min_y = [[sys.maxint for j in xrange(ydim)] for i in xrange(nBin)]
if err:
self.ySQlist = [[0.0 for j in xrange(ydim)] for i in xrange(nBin)]
def GetBin(self,x):
return int(floor((x-self.lo)/self.dx))
def whichBin(self,x,y, ibin=[] ):
assert x >= self.lo, 'x=%s,self.lo=%s'%(10**x,10**self.lo)
assert x < self.hi, 'x=%s,self.hi=%s'%(10**x,10**self.hi)
nr = int(floor((x-self.lo)/self.dx))
if ibin:
ibin[ 0 ] = nr
self.kounter[nr] += 1
self.xlist[nr] += x
for idim in xrange(len(y)):
self.ylist[nr][idim] += y[idim]
if y[idim] >= self.max_y[nr][idim]: #--- set max value
self.max_y[nr][idim]=y[idim]
if y[idim] <= self.min_y[nr][idim]:
self.min_y[nr][idim]=y[idim]
if self.err:
self.ySQlist[nr][idim] += y[ idim ] * y[ idim ]
def res(self, logScaleX = None, logScaleY = None, base = 10, SUM = None, MINMAX=None ):
indices = xrange(10**6)
someList = []
for x,index in zip(self.xlist,indices):
nb = self.kounter[index]
if nb == 0: continue
xbar = self.xlist[index]/nb
ybar = [y/nb for y in self.ylist[index]]
if self.err:
sigmaY = [ ( ysq / nb - YBAR * YBAR ) ** 0.5 / nb ** 0.5 for ysq, YBAR in zip( self.ySQlist[ index ], ybar )]
if SUM:
sigmaY = [ i * nb for i in sigmaY ]
if SUM:
ybar = [y for y in self.ylist[index]]
if MINMAX:
MAX_y = [y for y in self.max_y[index]]
if logScaleX:
xbar = base ** xbar
if logScaleY:
ybar = [ base ** item for item in ybar ]
if self.err:
someList.append([ xbar, ybar, sigmaY ])
elif MINMAX:
someList.append([ xbar, ybar, MAX_y ])
else:
someList.append([ xbar, ybar ])
return someList
class histogram( bins ):
def res( self, Radial = None, logScale = None, normalize = True, base = 10.0, ACCUMLTD = None ):
PDF = []
self.nPoint = nPoint = sum( self.kounter )
indices = xrange( sys.maxint )
y_accm = nPoint
for y, index in zip( self.kounter, indices ):
if not ACCUMLTD and y == 0:
continue
if not y == 0:
x = self.xlist[ index ] / y #self.lo + index * self.dx
else:
x = self.lo + index * self.dx
Y = 1.0 * y
dx = self.dx
if logScale:
x = base ** x
dx = x * ( base ** self.dx - 1.0 )
# print Y, dx
if normalize:
Y /= ( nPoint * dx )
if Radial:
Y /= ( 2.0 * pi * x )
# PDF.append( [ x, Y ] )
# PDF.append( [ x + dx, Y ] )
# PDF.append( [ x + 0.5 * dx, Y, 0.0, ( 1.0 * y_accm / nPoint if normalize else 1.0 * y_accm ) ] )
error_std = 0.0
if self.err:
error_std = sqrt( nPoint * Y * dx ) / ( nPoint * dx ) #--- poisson
error_std = sqrt( nPoint * Y * dx * ( 1.0 - Y * dx ) ) / ( nPoint * dx ) #--- bi-nomial
PDF.append( [ x, Y, 0.0, ( 1.0 * y_accm / nPoint if normalize else 1.0 * y_accm ), error_std ] )
y_accm -= y
return PDF
def GetCartesian( dff ):
df = dff.copy()
xlo = df['longitude'].min()
xhi = df['longitude'].max()
ylo = df['latitude'].min()
yhi = df['latitude'].max()
getDistX = lambda x: geopy.distance.vincenty( ( 0.0, xlo ), ( 0.0, x ) ).km
getDistY = lambda y: geopy.distance.vincenty( ( ylo, 0.0 ), ( y, 0.0 ) ).km
df[ 'x(km)' ] = df[ 'longitude' ].apply( getDistX )
df[ 'y(km)' ] = df[ 'latitude' ].apply( getDistY )
df[ 'z(km)' ] = df['depth']
return df
def fractalDimension2nd( coord ):
#--- sort
if type( coord ) == type( [] ):
coord = ListToDict( coord )
points = coord.keys()
points.sort()
hsObject = histogram( 18 * 8, log( 1e-10,10), log( 1e8, 10 ) )
for point_i in points:
for point_j in points:
if not point_i < point_j: #--- pairs ij with i<j
continue
rij = sum( [ ( i - j ) ** 2 for i, j in zip( coord[ point_i ], coord[ point_j ] ) ] ) # ** 0.5
assert rij > 0, 'rij=%s,coord[ %s ]=%s, coord[ %s ]=%s' %(rij,point_i,coord[ point_i ], point_j, coord[ point_j ] )
hsObject.whichBin( 0.5 * log( rij,10 ), [ 1.0 ] )
for items in hsObject.res( logScale = True, normalize = True, ACCUMLTD = True ):
if items[ 3 ] > 0.0:
yield items[ 0 ], items[ 3 ]
#---------------------------------------------------------------------------------
#-----------------
#-----------------
#-----------------
#---------------------------------------------------------------------------------
#--------------------
#----- subset
#--------------------
swarm_tmp = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,
myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_tmp, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi = swarm_lohi.sample( n = 2000 ) #--- sample
#--------------------
#--- cartesian coords
#--------------------
swarm_cartesian = GetCartesian( swarm_lohi )
#--------------------
#--- evaluate df
#--------------------
tmp_coord = swarm_cartesian['x(km)'].apply(lambda x: [x+rnd.random()*1e-6]) +\
swarm_cartesian['y(km)'].apply(lambda x: [x+rnd.random()*1e-6]) +\
swarm_cartesian['z(km)'].apply(lambda x: [x+rnd.random()*1e-6])
tmp_coord = tmp_coord.to_dict()
dict_NR = fractalDimension2nd( tmp_coord ) #, dmin = 1.0e-02 )
#--------------------
#--- scattered plot
#--------------------
swarm_cartesian.plot.scatter('longitude','latitude',
s=3**(swarm_lohi['magnitude']),
c='date',cmap='jet',
alpha = 0.4)
#--------------------
#--- N(r) vs r
#--------------------
plt.figure( figsize = (4,4))
plt.xlabel('r(km)')
plt.ylabel('N(r)')
plt.xlim(1e-3,1e5)
plt.ylim(1e-8,1)
plt.yscale('log')
plt.xscale('log')
d_f = 1.6
for i in dict_NR:
plt.plot([i[ 0 ]],
[1-i[ 1 ]],
'o',color='black') #--- accumulated
plt.plot(i[ 0 ],
(1-i[ 1 ])/i[0]**d_f,
'.',color='red') #--- accumulated
plt.show()
In [14]:
#--- evaluate fractal dimension (correltion dimension )
def getRmat2d( df_complete ):
nmax = len( df_complete )
rmat = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
for i in xrange( nmax ):
#--- distance matrix
df_dx = df_complete[ 'x(km)' ] - df_complete[ 'x(km)' ][ i ]
df_dy = df_complete[ 'y(km)' ] - df_complete[ 'y(km)' ][ i ]
df_sq = ( df_dx*df_dx+df_dy*df_dy ) ** 0.5
df_sq[ : i ] = 0
rmat[ i ] = np.array(df_sq)
return np.array( rmat )
def getRmat3d( df_complete ):
nmax = len( df_complete )
rmat = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
for i in xrange( nmax ):
#--- distance matrix
df_dx = df_complete[ 'x(km)' ] - df_complete[ 'x(km)' ][ i ]
df_dy = df_complete[ 'y(km)' ] - df_complete[ 'y(km)' ][ i ]
df_dz = df_complete[ 'z(km)' ] - df_complete[ 'z(km)' ][ i ]
df_sq = ( df_dx*df_dx+df_dy*df_dy+df_dz*df_dz ) ** 0.5
df_sq[ : i ] = 0
rmat[ i ] = np.array(df_sq)
return np.array( rmat )
#def GetCartesian( dff ):
# df = dff.copy()
# xlo = df['longitude'].min()
# xhi = df['longitude'].max()
# ylo = df['latitude'].min()
# yhi = df['latitude'].max()
# getDistX = lambda x: geopy.distance.vincenty( ( 0.0, xlo ), ( 0.0, x ) ).km
# getDistY = lambda y: geopy.distance.vincenty( ( ylo, 0.0 ), ( y, 0.0 ) ).km
# df[ 'r(km)' ] = df[ 'longitude' ].apply( getDistX ) + df[ 'latitude' ].apply( getDistY ) * 1j
# return df
# --------------------
# ----- subset
# --------------------
swarm_tmp = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,
myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_tmp, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
#--------------------
#--- cartesian coords
#--------------------
swarm_lohi = GetCartesian( swarm_lohi )
#swarm_lohi = DataFrameSubSet( swarm_lohi,
# 'x(km)',
# ( 180,300))
#swarm_lohi = DataFrameSubSet( swarm_lohi,
# 'y(km)',
# ( 100,250))
#--- we need the distance matrix!
#rmat2d = getRmat2d( swarm_lohi )
rmat3d = getRmat3d( swarm_lohi )
#r_list2d = rmat2d[rmat2d.nonzero()]
r_list3d = rmat3d[rmat3d.nonzero()]
# #--- histogram
r_min = log(np.min(r_list3d),10) #min(log(np.min(r_list2d),10),log(np.min(r_list3d),10))
r_max = log(np.max(r_list3d),10) #max(log(np.max(r_list2d),10),log(np.max(r_list3d),10))
nbins = int(ceil(r_max - r_min))*8
#2d
# hist2d, edges = np.histogram(r_list2d,
# bins=np.logspace(r_min, r_max, nbins) )
# r_bins, edges = np.histogram(r_list2d,
# bins=np.logspace(r_min, r_max, nbins), weights = r_list2d )
# r_bins2d = r_bins[hist2d != 0]
# hist2d = hist2d[hist2d != 0]
# r_bins2d /= hist2d
# # #--- normalize
# hist_acc_2d = np.cumsum(hist2d)*1.0
# hist_acc_2d /= hist_acc_2d[-1]
#3d
hist3d, edges = np.histogram(r_list3d,
bins=np.logspace(r_min, r_max, nbins) )
r_bins, edges = np.histogram(r_list3d,
bins=np.logspace(r_min, r_max, nbins), weights = r_list3d )
r_bins3d = r_bins[hist3d != 0]
hist3d = hist3d[hist3d != 0]
r_bins3d /= hist3d
#--- normalize
hist_acc_3d = np.cumsum(hist3d)*1.0
hist_acc_3d /= hist_acc_3d[-1]
##--------------------
##--- scattered plot
##--------------------
#plt.scatter(swarm_lohi['x(km)'],
# swarm_lohi['y(km)'],
# s=3**(swarm_lohi['magnitude']),
# cmap='jet',
# alpha = 0.1)
#
##--- plot
##d_f = 1.4
#fig = plt.figure(figsize=(4,4))
#ax = fig.add_subplot(111)
#ax.set_xlim(10**r_min,10**r_max)
#ax.set_xlabel('$r(km)$')
#ax.set_ylabel('$N(r)$')
#ax.set_xscale('log')
#ax.set_yscale('log')
#ax.errorbar( r_bins, hist_acc,
# yerr = hist_acc/hist**.5,
# color='black',marker='o')
#
#ax.errorbar( r_bins, 10000*hist_acc/r_bins**d_f, #--- rescaled
# yerr = hist_acc/hist**.5/r_bins**d_f,
# color='red',marker='o',label='$N/r^{%s}$'%d_f)
#ax.legend()
#
##fig.savefig('%s/corrDim.png'%DIR_OUTPT_figs,bbox_inches='tight')
In [15]:
#--- plot
d_f = 2.0
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.axis([1e-4,1e2,1e-6,1e0])
ax.set_xlabel('$r(km)$')
ax.set_ylabel('$N(r)$')
ax.set_xscale('log')
ax.set_yscale('log')
# ax.errorbar( r_bins2d, hist_acc_2d,
# yerr = hist_acc_2d/hist2d**.5,
# color='black',marker='o')
ax.errorbar( r_bins3d, hist_acc_3d,
hist_acc_3d/hist3d**.5,
color='black',marker='s',linestyle=None)
#ax.errorbar( r_bins, 10000*hist_acc/r_bins**d_f, #--- rescaled
# yerr = hist_acc/hist**.5/r_bins**d_f,
# color='red',marker='o',label='$N/r^{%s}$'%d_f)
ax.plot(r_bins3d, r_bins3d**d_f/10,
'-.r',label='$r^{%s}$'%d_f)
ax.legend()
fig.savefig('%s/corrDim.png'%DIR_OUTPT_figs,bbox_inches='tight')
plt.show()
In [16]:
#--- evaluate fractal dimension (box counting )
import geopy.distance
from math import *
import random as rnd
#import warnings
#warnings.filterwarnings('ignore') #--- get rid of warnings
#def GetCartesian( dff ):
# df = dff.copy()
# xlo = df['longitude'].min()
# xhi = df['longitude'].max()
# ylo = df['latitude'].min()
# yhi = df['latitude'].max()
# getDistX = lambda x: geopy.distance.vincenty( ( 0.0, xlo ), ( 0.0, x ) ).km
# getDistY = lambda y: geopy.distance.vincenty( ( ylo, 0.0 ), ( y, 0.0 ) ).km
# df[ 'r(km)' ] = df[ 'longitude' ].apply( getDistX ) + df[ 'latitude' ].apply( getDistY ) * 1j
# df[ 'x(km)' ] = df[ 'longitude' ].apply( getDistX )
# df[ 'y(km)' ] = df[ 'latitude' ].apply( getDistY )
# return df
#--------------------
#----- subset
#--------------------
# swarm_tmp = DataFrameSubSet( swarm,
# 'date',
# ( myTlohi.lo,
# myTlohi.hi ) )
# swarm_lohi = DataFrameSubSet( swarm_tmp, #--- complete catalog
# 'magnitude',
# ( mc, sys.maxint ) )
# #swarm_lohi = swarm_lohi.sample( n = 1000 ) #--- sample
# #--------------------
# #--- cartesian coords
# #--------------------
# swarm_lohi = GetCartesian( swarm_lohi )
# #swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
# # 'x(km)',
# # ( 200, 1000 ) )
# #swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
# # 'y(km)',
# # ( -1000, 220 ) )
# #--------------------
# #--- scattered plot
# #--------------------
# # plt.figure(figsize=(8,8))
# # plt.xlim(0,60)
# # plt.ylim(-15,30+15)
# # plt.xlabel('x(km)')
# # plt.ylabel('z(km)')
# swarm_cp = swarm_lohi[swarm_lohi['magnitude']>2]
# swarm_cp.plot.scatter('x(km)','z(km)',
# s=3**(swarm_cp['magnitude']),
# cmap='jet',
# c='magnitude' )
# # alpha = 0.01)
# #--------------------
# #--- box counting
# #--------------------
# #tmp=pd.DataFrame({'x(m)':swarm_lohi['y(m)'],'y(m)':swarm_lohi['z(m)']}) #,'z(m)':swarm_lohi['z(m)']})
# #rmat = getRmat2d( tmp )
# #rmat = getRmat2d( swarm_lohi )
# #r_list = rmat[rmat.nonzero()]
# #r_min = np.min(r_list)
# #r_max = np.max(r_list)
# #
# #plt.savefig('%s/map2d.png'%DIR_OUTPT_figs,bbox_inches='tight')
In [17]:
#--- evaluate fractal dimension (box counting )
def fractalDimension3d( x, y, z, L, dmin = None ):
x = np.array(x)
y=np.array(y)
z=np.array(z)
# xlo = x.min()
# xhi = x.max()
# ylo = y.min()
# yhi = y.max()
# zlo = z.min()
# zhi = z.max()
# L = max( zhi - zlo, yhi - ylo, xhi - xlo )
#---
xc = x.mean() #--- center
yc = y.mean()
zc = z.mean()
xo = xc - L #--- origin
yo = yc - L
zo = zc - L
# L *= 2.0
#---
kounter = 0 #--- number of divisions
dx = dy = dz = L #--- scale
dl = []
N = []
while dx - dmin > 0.0: #--- do until dx < dmin
d = {}
index_i = np.array( map(int,np.floor( ( y - yo ) / dy ) ) )
index_j = np.array( map( int, np.floor( ( x - xo ) / dx ) ) )
index_k = np.array( map( int, np.floor( ( z - zo ) / dz ) ) )
inside_cell_lo = np.all([index_i>=0,index_j>=0,index_k>=0.0],axis=0) #--- point is inside?
inside_cell_hi = np.all([index_i<2**kounter,index_j<2**kounter,index_k<2**kounter],axis=0) #
inside_cell = np.all([inside_cell_lo,inside_cell_hi],axis = 0)
cell_id = index_k * 4 ** kounter + index_i * ( 2**kounter ) + index_j
#---
dl+=[dx]
N+=[len( set( cell_id ) )]
#---
kounter += 1
dx = dy = dz = L / 2**kounter
assert dx - dmin <= 0.0
#---
return np.array(dl), np.array(N)
def fractalDimension2d( x, y, L, dmin = None ):
x = np.array(x)
y=np.array(y)
# xlo = x.min()
# xhi = x.max()
# ylo = y.min()
# yhi = y.max()
# L = max( yhi - ylo, xhi - xlo )
#---
xc = x.mean() #--- center
yc = y.mean()
xo = xc - L #--- origin
yo = yc - L
# L *= 2.0
#---
kounter = 0 #--- number of divisions
dx = dy = L / 2**kounter #--- scale
dl = []
N = []
# pdb.set_trace()t
# fig = plt.figure(figsize=(4,4))
# ax=fig.add_subplot(111)
# plt.scatter(x,y)
# ax.add_patch( patches.Rectangle(xy=(xo,yo), width=L,
# height=L, linewidth=1,
# clip_on=False,facecolor=None,edgecolor='black',fill=None) )
while dx - dmin > 0.0: #--- do until dx < dmin
d = {}
index_i = np.array( map(int,np.floor( ( y - yo ) / dy ) ) ) #--- assingn index
index_j = np.array( map( int, np.floor( ( x - xo ) / dx ) ) )
inside_cell_lo = np.all([index_i>=0,index_j>=0],axis=0) #--- point is inside?
inside_cell_hi = np.all([index_i<2**kounter,index_j<2**kounter],axis=0) #
inside_cell = np.all([inside_cell_lo,inside_cell_hi],axis = 0)
index_i = index_i[ inside_cell ] #--- subset of the index array
index_j = index_j[ inside_cell ]
cell_id = index_i * ( 2**kounter ) + index_j
#---
dl+=[dx]
N+=[len( set( cell_id ) )]
#---
kounter += 1
dx = dy = L / 2**kounter
assert dx - dmin <= 0.0
#---
return np.array(dl), np.array(N)
# fig = plt.figure(figsize=(10,4))
# ax = fig.add_subplot(121)
# ax2 = fig.add_subplot(122)
# ax.set_xlabel('$r(km)$')
# ax.set_ylabel('$N(r)$')
# ax.set_xscale('log')
# ax.set_yscale('log')
# ax2.set_xlabel('$r(km)$')
# ax2.set_ylabel(r'$N(r)\times r^{%s}$'%d_f)
# ax2.set_xscale('log')
# ax2.set_yscale('log')
# dl={}
# N={}
# count = 0
# for L in np.logspace(log(r_min,2),log(r_max,2),num=6,base=2.0):
# #for L in [2*420.0]:
# dl[count], N[count] = fractalDimension2d( swarm_lohi['x(km)'],
# swarm_lohi['y(km)'], L, dmin=r_min ) # , swarm_lohi['z(m)'], dmin=r_min )
# # dl[count], N[count] = fractalDimension3d( swarm_lohi['x(m)'],
# # swarm_lohi['y(m)'],
# # swarm_lohi['z(m)'],
# # L, dmin=r_min )
# ##--- try random
# #x=np.array([rnd.random() for i in xrange(1000)])
# #y=0.2*x #[rnd.random() for i in xrange(10000)]
# ##z=x #[rnd.random() for i in xrange(10000)]
# #tmp=pd.DataFrame({'x(m)':x,'y(m)':y})#,'z(m)':z})
# #rmat = getRmat2d( tmp )
# #r_list = rmat[rmat.nonzero()]
# #r_min = min(r_list)
# #sdict = fractalDimension2d( x,y,
# # dmin=r_min )
# #plt.scatter(x,y)
# #--- plot
# #ax.set_xlim(10**r_min,10**r_max)
# #for key in sdict:
# ax.plot(dl[count],N[count],'-o')
# ax2.plot(dl[count],N[count]*dl[count]**d_f,marker='o',label='$L=%2.1e$'%L)#,linestyle='None')
# count+=1
# ax2.legend(bbox_to_anchor=(1, 1))
# fig.tight_layout()
# #fig.savefig('%s/boxCount.png'%DIR_OUTPT_figs,bbox_inches='tight')
In [18]:
# comment if passed by arg
#BVALL = 0.84
#MCC = 1.5
Df = 2.0
bval = BVALL
mc = MCC
#print BVAL, MC
#print bval, mc
In [19]:
# --- trig analysis (new) (it's a large data sets. run the simulation on a cluster and use the outputs)
from scipy.sparse import lil_matrix
import time
from IPython.display import display
import datetime
import pdb
import numpy.linalg
import numpy as np
def UpperTriang(a):
il1 = np.tril_indices(a.shape[0])
a[il1] = 0
return a
def getTmat( df_complete ):
nmax = len( df_complete )
prefact = 1.0 / ( 24.0 * 60 * 60 ) #--- daily
tmat = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
tmat_ut = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
for i in xrange( nmax ):
df_dt = ( df_complete[ 'date' ] - df_complete[ 'date' ][ i ] ).dt.total_seconds() * prefact #--- time diff between i-th event and all subsequent events
tmat[ i ] = np.array(df_dt)
df_dt[ : i ] = df_dt[ i ] #--- must have an upper triangular matrix
tmat_ut[ i ] = np.array(df_dt)
#---
return np.array( tmat ), np.array( tmat_ut )
def AddPositionList(df_complete):
df_complete['r(m)'] = df_complete['x(m)'].apply(lambda x: [x])+ \
df_complete['y(m)'].apply(lambda x: [x])+ \
df_complete['z(m)'].apply(lambda x: [x])
def getRmat( df_complete ):
nmax = len( df_complete )
rmat = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
rmat_ut = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
for i in xrange( nmax ):
#--- distance matrix
df_dx = df_complete[ 'x(km)' ] - df_complete[ 'x(km)' ][ i ]
df_dy = df_complete[ 'y(km)' ] - df_complete[ 'y(km)' ][ i ]
df_dz = df_complete[ 'z(km)' ] - df_complete[ 'z(km)' ][ i ]
df_sq = ( df_dx*df_dx+df_dy*df_dy+df_dz*df_dz ) ** 0.5
rmat[ i ] = np.array(df_sq)
df_sq[ : i ] = 0
rmat_ut[ i ] = np.array(df_sq)
return np.array( rmat ), np.array( rmat_ut )
def getMmat( df_complete ):
nmax = len( df_complete )
m_mat = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
m_mat_ut = np.matrix(np.zeros(nmax*nmax).reshape(nmax,nmax))
for i in xrange( nmax ):
#--- magnitude
df_m = pd.Series( [ df_complete[ 'magnitude' ][ i ] ] * nmax )
m_mat[ i ] = np.array(df_m)
df_m[ : i ] = 0
m_mat_ut[ i ] = np.array(df_m)
return np.array( m_mat ), np.array( m_mat_ut )
def vectorizedAnalysis( df_complete ):
#--- setup t, r, m
t0 = time.time()
tmat0, tmat = getTmat( df_complete )
print 'setting up tmat:%s s'%(time.time() - t0)
t0 = time.time()
r_mat0, r_mat = getRmat( df_complete )
print 'setting up rmat:%s s'%(time.time() - t0)
t0 = time.time()
m_mat0, m_mat = getMmat( df_complete )
print 'setting up m_mat:%s s'%(time.time() - t0)
# #--- make upper-triangular matrices
# t0 = time.time()
# tmat = UpperTriang(np.copy(tmat0))
# r_mat = UpperTriang(np.copy(r_mat0))
# m_mat = UpperTriang(np.copy(m_mat0))
# print 'setting up ut matrix:%s s'%(time.time() - t0)
#--- nij
t0 = time.time()
NIJ = tmat * r_mat ** ( Df ) * 10 ** ( - bval * ( m_mat - mc ) )
TIJ = tmat * 10 ** ( - 0.5 * bval * ( m_mat - mc ) ) #--- scaled time
RIJ = r_mat ** ( Df ) * 10 ** ( - 0.5 * bval * ( m_mat - mc ) ) #--- scaled time
print 'setting up NIJ:%s s'%(time.time() - t0)
nmax = len( df_complete )
N_sparse = lil_matrix( ( nmax, nmax ) ) #--- sparse matrix
T_sparse = lil_matrix( ( nmax, nmax ) )
R_sparse = lil_matrix( ( nmax, nmax ) )
for junk, j in zip( NIJ, xrange( sys.maxint ) ):
if j == 0:
continue
x = min( NIJ[ : j, j ][ NIJ[ : j, j ] != 0 ] ) #--- min nij
assert x != 0.0
rowi = np.where( NIJ[ : j, j ] == x )[ 0 ][ 0 ] #--- find row
assert rowi <= j
N_sparse[ rowi, j ] = x #--- insert
T_sparse[ rowi, j ] = TIJ[ rowi, j ]
R_sparse[ rowi, j ] = RIJ[ rowi, j ]
return N_sparse, (T_sparse, tmat0 ), ( R_sparse, r_mat0 ), m_mat0
#-------------------
#--- set parameters
#-------------------
# Df = 2.0
# bval = 0.95
# mc=1.0
#--------------------
#----- subset
#--------------------
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi.reset_index(inplace=True,drop=True)
#--------------------
#--- cartesian coords
#--------------------
swarm_lohi = GetCartesian( swarm_lohi )
#AddPositionList( swarm_lohi )
swarm_lohi.reset_index(inplace=True,drop=True)
#--- vectorized
t0 = time.time()
N_sparse, (T_sparse,tmat), (R_sparse,rmat), m_mat = vectorizedAnalysis( swarm_lohi )
print 'duration=%s s'%(time.time() - t0)
In [20]:
#--- random analysis
def shuffle(amat):
ind_list = range(amat.shape[0])
np.random.shuffle(ind_list)
amat_row = amat[ind_list]
return amat_row[:,ind_list]
def GetNij(rmat,tmat,m_mat,df_complete):
tmat = shuffle(np.abs(tmat))
r_mat = shuffle(rmat)
m_mat = shuffle(m_mat)
#--- nij
NIJ = tmat * r_mat ** ( Df ) * 10 ** ( - bval * ( m_mat - mc ) )
TIJ = tmat * 10 ** ( - 0.5 * bval * ( m_mat - mc ) ) #--- scaled time
RIJ = r_mat ** ( Df ) * 10 ** ( - 0.5 * bval * ( m_mat - mc ) ) #--- scaled time
nmax = len( df_complete )
N_sparse = lil_matrix( ( nmax, nmax ) ) #--- sparse matrix
T_sparse = lil_matrix( ( nmax, nmax ) )
R_sparse = lil_matrix( ( nmax, nmax ) )
for junk, j in zip( NIJ, xrange( sys.maxint ) ):
if j == 0:
continue
x = min( NIJ[ : j, j ][ NIJ[ : j, j ] != 0 ] ) #--- min nij
assert x != 0.0
rowi = np.where( NIJ[ : j, j ] == x )[ 0 ][ 0 ] #--- find row
N_sparse[ rowi, j ] = x #--- insert
T_sparse[ rowi, j ] = TIJ[ rowi, j ]
R_sparse[ rowi, j ] = RIJ[ rowi, j ]
return N_sparse, (T_sparse, tmat ), ( R_sparse, r_mat), m_mat
N_sparse_rnd = {}
T_sparse_rnd = {}
R_sparse_rnd = {}
tmat_rnd = {}
rmat_rnd = {}
nsamp = 1
tmat = np.abs(tmat)
if 1:
#for i in xrange(nsamp):
# print i
N_sparse_rnd, (T_sparse_rnd,tmat_rnd), (R_sparse_rnd,rmat_rnd), junk = GetNij(rmat,tmat,m_mat,swarm_lohi)
In [21]:
#--- save sparse matrices
import scipy.sparse
for mat, title in zip([N_sparse,T_sparse,R_sparse],
['N_sparse_matrix.npz','T_sparse_matrix.npz','R_sparse_matrix']):
scipy.sparse.save_npz('%s/%s'%(DIR_OUTPT,title), mat.tocsr()) #--- output
for mat, title in zip([N_sparse_rnd,T_sparse_rnd,R_sparse_rnd],
['N_sparse_rnd_matrix.npz','T_sparse_rnd_matrix.npz','R_sparse_rnd_matrix']):
scipy.sparse.save_npz('%s/%s'%(DIR_OUTPT,title), mat.tocsr()) #--- output
#swarm_shuffled.to_csv('%s/swarm_shuffled.csv'%DIR_OUTPT)
In [22]:
#--- load sparse matrices
#--- 1- run on the cluster:
#--- 2- sbatch --mem=8gb --partition=single --time=02:59:59 -n 1 ./oarScript.sh
#--- 3- oarScript.sh: source activate conda-env; python ./swarmEMC.py
#--- 4- copy: scp arc:/home/kamran.karimi1/Project/seismic/dataset/El_Mayor_Cucpah/* ./swarm/dataset/El_Mayor_Cucpah/
import scipy.sparse
N_sparse = scipy.sparse.load_npz('%s/N_sparse_matrix.npz'%DIR_OUTPT )
T_sparse = scipy.sparse.load_npz('%s/T_sparse_matrix.npz'%DIR_OUTPT )
R_sparse = scipy.sparse.load_npz('%s/R_sparse_matrix.npz'%DIR_OUTPT )
#rmat = np.load('%s/rmat.npy'%DIR_OUTPT)
#tmat = np.load('%s/tmat.npy'%DIR_OUTPT)
nsamp = 3
#N_sparse_rnd = {}
#T_sparse_rnd = {}
#R_sparse_rnd = {}
if 1:
#for i in xrange(nsamp):
N_sparse_rnd = scipy.sparse.load_npz('%s/N_sparse_rnd_matrix.npz'%(DIR_OUTPT ) )
T_sparse_rnd = scipy.sparse.load_npz('%s/T_sparse_rnd_matrix.npz'%(DIR_OUTPT ) )
R_sparse_rnd = scipy.sparse.load_npz('%s/R_sparse_rnd_matrix.npz'%(DIR_OUTPT ) )
#rmat_rnd = np.load('%s/rmat_rnd.npy'%DIR_OUTPT)
#tmat_rnd = np.load('%s/tmat_rnd.npy'%DIR_OUTPT)
#swarm_shuffled = pd.read_csv('%s/swarm_shuffled.csv'%DIR_OUTPT)
#swarm_shuffled_incomplete = pd.read_csv('%s/swarm_shuffled_incomplete.csv'%DIR_OUTPT)
In [23]:
##--- matrix map of nij & nij_rand
#
#plt.subplot(2,1,1)
#plt.pcolormesh(np.log(N_sparse.toarray()),cmap='jet')
#plt.subplot(2,1,2)
#plt.pcolormesh(np.log(N_sparse_rnd.toarray()),cmap='jet')
In [24]:
mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings
rcc('text', usetex=True)
font = {'size' : 26}
matplotlib.rc('font', **font)
#---------------------------
#--- compute threshold n_th
#---------------------------
def PLOT(edges, hist, cdf, title ):
#--- plot
plt.figure(figsize=(8,4))
plt.subplot(1,2,1)
# plt.xlim(1e-13,1)
# plt.ylim(1e-1,1e13)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('nij')
plt.ylabel('P(nij)')
plt.title( title )
plt.plot(edges[:-1],hist,'.-')
plt.subplot(1,2,2)
# plt.xlim(1e-13,1)
# plt.ylim(1e-3,1)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('nij')
plt.ylabel('cdf')
plt.plot(edges[:-1],cdf,'.-')
plt.plot([n_thresh for i in edges[:-1]], cdf)
#---
def HIST( nij, nbins ):
#--- histograms
nij_TxR = nij['T*']*nij['R*']
hist, edges = np.histogram( nij_TxR,
bins=np.logspace(np.log(nij_TxR.min()),np.log(nij_TxR.max()),nbins),
density=True)
#--- accumulated
cdf, edges = np.histogram( nij_TxR,
bins=np.logspace(np.log(nij_TxR.min()),np.log(nij_TxR.max()),nbins,endpoint=True))
cdf=np.cumsum(cdf)
#assert cdf[-1] == len(nij_TxR), '%s,%s'%(cdf[-1], len(nij_TxR))
N = 1.0 * len(nij_TxR)
cdf = cdf / N
return edges, hist, cdf
def reset_xticks(nn,ax,(TMIN,TMAX)):
ax.set_xticks( np.linspace(ax.axis()[0],ax.axis()[1],nn) )
labels=np.linspace(TMIN,TMAX,nn)
labels = [r'$10^{%i}$'%i for i in labels]
ax.set_xticklabels(labels)
def reset_yticks(nn,ax,(RMIN,RMAX)):
ax.set_yticks( np.linspace(ax.axis()[2],ax.axis()[3],nn) )
labels=np.linspace(RMIN,RMAX,nn)
labels = [r'$10^{%i}$'%i for i in labels]
ax.set_yticklabels(labels)
def GetTrigListLoHi( nij_trig, col, (tlo,thi), swarm_lohi ):
list_of_mothers = nij_trig.groupby(by='Parent_id').groups.keys()
list_of_mothers.sort()
tmp_df = swarm_lohi.loc[ list_of_mothers ]
tmp_df = tmp_df[ (tmp_df[ col ] >= tlo) &
(tmp_df[ col ] < thi) ]
return nij_trig[ pd.DataFrame([nij_trig['Parent_id'] == i for i in tmp_df.index]).any() ]
def GetInterpolatedData( df0, (TMIN,TMAX),(RMIN,RMAX) ):
nbins_per_decade = 16
nbins_x=nbins_y=int(TMAX-TMIN)*nbins_per_decade
df=df0.copy()
df['T*']=df[ 'T*' ].apply(lambda x: log(x,10))
df['R*']=df[ 'R*' ].apply(lambda x: log(x,10))
heatmap, xedges, yedges = np.histogram2d( df[ 'T*' ], df[ 'R*'], bins=[np.linspace(TMIN,TMAX,nbins_x+1), np.linspace(RMIN,RMAX,nbins_y+1)], normed=True)
heatmap *= len( df )
heatmap = gaussian_filter( heatmap, sigma = nbins_per_decade/4 )
return heatmap
def GetDistance( df1, df2 ):
df_sq = (np.array( df1['x(km)'] ) - np.array( df2['x(km)'] ) ) ** 2 + \
(np.array( df1['y(km)'] ) - np.array( df2['y(km)'] ) ) ** 2 + \
(np.array( df1['z(km)'] ) - np.array( df2['z(km)'] ) ) ** 2
return pd.Series( df_sq ** 0.5 )
def RemovePair( dff, cs, catalog ):
pref = 1.0 / 3600.0 / 24 #--- daily
df = dff.copy()
series_r = GetDistance(catalog.loc[df['Parent_id']],catalog.loc[df['Event_id']])
# series_r = pd.Series(np.array( catalog.loc[df['Parent_id']]['r(km)'] ) -\
# np.array( catalog.loc[df['Event_id']]['r(km)'] ) )
assert len(series_r[series_r==np.nan]) == 0
df['R'] = series_r.copy() #abs()
series_t = pd.Series(np.array( catalog.loc[df['Event_id']]['date'] ) -\
np.array( catalog.loc[df['Parent_id']]['date'] ) )
df['T'] = series_t.apply(lambda x:x.total_seconds()*pref)
assert len ( df[ df['R'] == 0.0 ] ) == 0, '%s'%display( df[ df['R'] == 0.0 ] )
assert len ( df[ df['T'] == 0.0 ] ) == 0, '%s'%display( df[ df['T'] == 0.0 ] )
return df [ df[ 'R' ] <= df[ 'T' ] * cs ]
#--------------------
#----- set parameters
#--------------------
#--- set quantile
#quantile = 0.05 #new_quantile #0.050
#--- wave speed
cs = 3 * (24*3600) #--- km per day
#--- cut-off time
dt_cut = 0.0 #--- day
#--- plot scattered
TMIN=-5
TMAX=3
RMIN=-4
RMAX=4
#--------------------
#----- subset
#--------------------
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi.reset_index(inplace=True,drop=True)
#--------------------
#--- cartesian coords
#--------------------
swarm_lohi = GetCartesian( swarm_lohi )
#AddPositionList( swarm_lohi )
swarm_lohi.reset_index(inplace=True,drop=True)
#--------------------
#--- store actual data in a dataframe
#--------------------
rows, cols = N_sparse.nonzero()
tlist = T_sparse[rows, cols].tolist()[0]
rlist = R_sparse[rows, cols].tolist()[0]
nij = pd.DataFrame({'Event_id':cols, 'Parent_id':rows,'T*':tlist,'R*':rlist})
#--------------------
#--- store random data in a dataframe
#--------------------
#nij_rnd = {}
#nsamp = 1
#for i in xrange(nsamp):
i=0
rows, cols = N_sparse_rnd.nonzero()
tlist = T_sparse_rnd[rows, cols].tolist()[0]
rlist = R_sparse_rnd[rows, cols].tolist()[0]
nnn = len(cols)
nij_rnd = pd.DataFrame({'Event_id':cols, 'Parent_id':rows,'T*':tlist,'R*':rlist},index=range(i*nnn,(i+1)*nnn))
#--------------------------------------------
#----- only include mother events with
#----- t_mother between (t_lo, t_hi)
#--------------------------------------------
#--- actual data
#--- remove pairs with r > cs.t
nij = RemovePair( nij, cs, swarm_lohi )
myTlohi2nd = Tlohi(myTlohi.lo+datetime.timedelta(days=dt_cut), #--- set cut-off time
myTlohi.hi+datetime.timedelta(days=-dt_cut)
)
nij = GetTrigListLoHi( nij, 'date', #--- t_mother > t_lo
( myTlohi2nd.lo, myTlohi2nd.hi ),
swarm_lohi)
nij.reset_index(inplace=True, drop=True)
#--- random data
#for i in xrange(nsamp):
# nij_rnd = RemovePair( nij_rnd, cs, swarm_shuffled )
# nij_rnd = GetTrigListLoHi( nij_rnd, 'date',
# ( myTlohi.lo, myTlohi.hi ),
# swarm_shuffled)
# nij_rnd.reset_index(inplace=True, drop=True)
#--- compute threshold
N_thresh = []
frames = []
#for i in xrange(nsamp):
# frames.append(nij_rnd)
#--- plot histograms for random data
edges, hist, cdf = HIST( nij_rnd, 1024 )
# edges = 0.5*(edges[:-1]+edges[1:])
#n_thresh = edges[:-1][cdf>quantile][0]
#n_thresh = np.mean(N_thresh)
#print 'n_th=%2.1e+%2.1e' %(n_thresh,np.std(N_thresh) / nsamp ** 0.5)
Nij_rnd = nij_rnd #pd.concat(frames)
#--- distribution of nij (actual data)
edges_actual, hist_actual, cdf_actual = HIST( nij, 1024 )
#edges_actual = 0.5*(edges_actual[:-1]+edges_actual[1:])
#PLOT(edges, hist, cdf, 'random')
Evaluate threshold: Gaussian Mixture Model¶
In [25]:
import sklearn
from sklearn import linear_model, mixture
import sklearn.mixture.gaussian_mixture as gm
import sklearn.mixture.bayesian_mixture as bm
from scipy import optimize
import scipy
def gaussian(x, mu, sig):
return 1./(np.sqrt(2.*np.pi)*sig)*np.exp(-np.power((x - mu)/sig, 2.)/2)
def gaussian_mixture( nij, edge_act, hist_act, LABELS = False):#,
# edge_red, hist_red,
# edge_blue, hist_blue, nmin, nmax):
#--- convert data frema to array
X=np.log10(nij['R*']*nij['T*'])
X=np.array(X).reshape(len(X),1)
#--- initilize classifier
# gm_obj = gm.GaussianMixture(n_components=2, tol=1e-14, max_iter=1000000,verbose=1)
# gm_obj.fit(X)
gm_obj = bm.BayesianGaussianMixture(n_components=2, tol=1e-08, max_iter=1000000,
verbose=0,reg_covar=1e-08,weight_concentration_prior=10000.0) #,weight_concentration_prior_type = 'dirichlet_distribution')
gm_obj.fit(X)
#--- min(\mu0, \mu1) corresponds to trig. mode
mean0 = gm_obj.means_[0][0]
mean1 = gm_obj.means_[1][0]
d = { 'trigrd':min([mean0,0],[mean1,1])[1],
'backgrd':max([mean0,0],[mean1,1])[1]}
mean = { 'trigrd' : gm_obj.means_[d['trigrd']][0],
'backgrd': gm_obj.means_[d['backgrd']][0] }
sigma = { 'trigrd' : gm_obj.covariances_[d['trigrd']][0]**0.5,
'backgrd': gm_obj.covariances_[d['backgrd']][0]**0.5 }
delta = { 'trigrd' : gm_obj.weights_[d['trigrd']],
'backgrd': gm_obj.weights_[d['backgrd']] }
#--- plot scatter
nij_red = nij[gm_obj.predict(X)==0]
nij_blue = nij[gm_obj.predict(X)==1]
prob = np.random.rand(len(X))
list_of_red = prob < gm_obj.predict_proba(X)[:, d['trigrd']] #--- assign modes based on the responsibilities \gamma
list_of_blue = prob >= gm_obj.predict_proba(X)[:, d['trigrd']]
nij_red = nij[list_of_red]
nij_blue = nij[list_of_blue]
#
plotScatter(nij_red, nij_blue,
# TMIN=-13.0,TMAX=-4.0,RMIN=-6.0,RMAX=3.0,
label0='Trig.',label1='Backgrd',
save = '%s/densityMapsScatterGM.png'%DIR_OUTPT_figs )
#--- plot distributions
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.set_yscale('log')
ax.set_ylim(1e0,1e3)
#
ax.plot(edge_act,hist_act,'o',color='black',label='Total')
# ax.plot(edge_red,hist_red,'s',color='red',label='trig',alpha=0.4)
# ax.plot(edge_blue, hist_blue,'^',color='C0',label='back',alpha=0.4)
# #
xv = edge_act #np.linspace(np.min(X),np.max(X),128)
ax.plot( xv,
len(X)*(delta['trigrd']*gaussian(xv, mean['trigrd'], sigma['trigrd'])+
delta['backgrd']*gaussian(xv, mean['backgrd'], sigma['backgrd']))*(xv[1]-xv[0]),
color='black')
ax.plot( xv, len(X)*delta['trigrd']*gaussian(xv, mean['trigrd'], sigma['trigrd'])*(xv[1]-xv[0]),color='red')
ax.plot( xv, len(X)*delta['backgrd']*gaussian(xv, mean['backgrd'], sigma['backgrd'])*(xv[1]-xv[0]),color='C0')
if LABELS:
ax.set_xlabel(r'$log n_{ij}$')
ax.set_ylabel(r'$P(log n_{ij})$')
ax.set_xlim(floor(np.min(edge_act)),ceil(np.max(edge_act)))
#
fig.savefig('%s/distnijGM.png'%DIR_OUTPT_figs,bbox_inches='tight',dpi=2*75)
plt.show()
#--- find decision boundary
mu0,mu1,sigma0,sigma1,delta0,delta1=mean['trigrd'],mean['backgrd'],sigma['trigrd'],sigma['backgrd'],delta['trigrd'],delta['backgrd']
def f(x): #,:
return delta0*np.exp(-0.5*((x-mu0)/sigma0)**2)/sigma0 - delta1*np.exp(-0.5*((x-mu1)/sigma1)**2)/sigma1
n_th = optimize.bisect(f, mu0-1*sigma0, mu1+1*sigma1)
return n_th
def DistNij(nij,normed=True, nbins_per_decade = 4, **kwargs):
#--- histogram
nn=np.log10(nij['R*']*nij['T*'])
nmin=kwargs['nmin'] if 'nmin' in kwargs else nn.min()
nmax=kwargs['nmax'] if 'nmax' in kwargs else nn.max()
bins=np.linspace(nmin,nmax,int(nmax-nmin)*nbins_per_decade)
hist, edge = np.histogram(nn,bins=bins,normed=normed)
#--- accumulated histogram
slist=np.array(nn)
slist.sort()
N = len(slist)
d = histogramACCUMLTD( slist.tolist() )
keys=d.keys()
keys.sort()
xx=[];yy=[]
for ikey in keys:
xx.append(d[ikey][0])
yy.append(d[ikey][2])
# ax2 = ax.twinx()
# ax2.plot(xx,yy,
# linestyle='-', drawstyle='steps-post',color='red',
# linewidth=1.0) #--- accumulated
# #ax2.set_xlim(-7,1)
# #ax2.set_ylim(0,1200)
# ax2.tick_params(axis='y',colors='red')
# ax2.set_ylabel('$N(<n_{ij})$',color='red')
return edge[:-1],hist
def plotScatter(nij, nij_rnd, **kwargs):
#--- (Tmin,Tmax)-(Rmin,Rmax)
TMIN = kwargs['TMIN'] if 'TMIN' in kwargs else np.floor(np.log10(nij['T*']).min())
TMAX = kwargs['TMAX'] if 'TMAX' in kwargs else np.ceil(np.log10(nij['T*']).max())
RMIN = kwargs['RMIN'] if 'RMIN' in kwargs else np.floor(np.log10(nij['R*']).min())
RMAX= kwargs['RMAX'] if 'RMAX' in kwargs else np.ceil(np.log10(nij['R*']).max())
label=['','']
label[0] = kwargs['label0'] if 'label0' in kwargs else ''
label[1] = kwargs['label1'] if 'label1' in kwargs else ''
color=['red','C0']
fig = plt.figure(figsize=(4,4))
for nn, index, title in zip([ nij, nij_rnd ],xrange(2),['actual','randomized']):
ax = fig.add_subplot(111) #,2,index+1)
#--- interpolate
reset_xticks(4,ax,(TMIN,TMAX))
#--- reset y-ticks
reset_yticks(4,ax,(RMIN,RMAX))
#--- nij must be within limits
nij_copy = nn.copy()
nij_copy=nij_copy[(nij_copy['T*'] >= 10**TMIN ) &
(nij_copy['T*'] < 10**TMAX ) &
(nij_copy['R*'] >= 10**RMIN ) &
(nij_copy['R*'] < 10**RMAX )]
#--- mapping
xx=ax.axis()[0]+(ax.axis()[1]-ax.axis()[0])*(np.log10(nij_copy['T*'])-TMIN)/(TMAX-TMIN)
yy=ax.axis()[2]+(ax.axis()[3]-ax.axis()[2])*(np.log10(nij_copy['R*'])-RMIN)/(RMAX-RMIN)
ax.set_xlim(ax.axis()[0],ax.axis()[1])
ax.set_ylim(ax.axis()[2],ax.axis()[3])
#--- plot
# DrawFrame(ax, (0.25,0.08),(0.16,0.08),0.01)
ax.scatter(xx,
yy,
alpha=1,color=color[index],s=0.1,label=label[index])
ax.legend(fontsize='xx-small')
# #--- draw the separation line
# xx=ax.axis()[0]+(ax.axis()[1]-ax.axis()[0])*(np.log10(nij_copy['T*'])-TMIN)/(TMAX-TMIN)
# yy=ax.axis()[2]+(ax.axis()[3]-ax.axis()[2])*(np.log10(n_thresh/nij_copy['T*'])-RMIN)/(RMAX-RMIN)
ax.set_xlabel(r'$T$')
ax.set_ylabel(r'$R$')
# ax.plot( xx, yy,'-.',color='white',linewidth=1)
if 'save' in kwargs:
fig.savefig(kwargs['save'],bbox_inches='tight',dpi=75)
plt.show()
#--- actual
edge_act, hist_act = DistNij(nij,normed=None, nbins_per_decade = 8)
n_thresh = gaussian_mixture( nij, edge_act, hist_act)#,
# edge_red, hist_red,
# edge_blue, hist_blue, nmin, nmax)
n_thresh = 10 ** n_thresh
Distributions of $n_{ij}$¶
In [26]:
#--- histogram
nn=np.log10(nij['R*']*nij['T*'])
nmin=nn.min()
nmax=nn.max()
bins=np.linspace(nmin,nmax,int(nmax-nmin)*6)
hist, edge = np.histogram(nn,bins=bins,normed=True)
count, junk = np.histogram(nn,bins=bins)
hist=hist[count!=0]
edge=edge[:-1][count!=0]
count = count[ count !=0 ]
err_hist = hist / count**.5
#--- accumulated histogram
slist=np.array(nn)
slist.sort()
N = len(slist)
d = histogramACCUMLTD( slist.tolist() )
keys=d.keys()
keys.sort()
xx=[];yy=[]
for ikey in keys:
xx.append(d[ikey][0])
yy.append(d[ikey][2])
#------------
#--- plot
#------------
fig= plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.set_xlabel('$log~n_{ij}$')
ax.set_ylabel('$P(log~n_{ij})$')
(nmin,nmax) = (min(slist),max(slist))
ax.set_xlim(floor(nmin),ceil(nmax))
ax.set_xlim(-7,1)
ax.set_ylim(1e-2,1)
#ax.set_xscale('log')
ax.set_yscale('log')
ax.errorbar(edge,hist,yerr=err_hist,color='black',
fmt='-o', markersize=8,markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1)
# ax2 = ax.twinx()
# ax2.plot(xx,yy,
# linestyle='-', drawstyle='steps-post',color='red',
# linewidth=1.0) #--- accumulated
# ax2.set_xlim(-7,1)
# ax2.set_ylim(0,1200)
# ax2.tick_params(axis='y',colors='red')
# ax2.set_ylabel('$N(<n_{ij})$',color='red')
np.save('actual',np.c_[edge,hist,err_hist])
#np.save('acc_actual',np.c_[xx,yy])
#fig.savefig('%s/pn_actual.png'%DIR_OUTPT_figs,bbox_inches='tight')
plt.show()
Distributions of random $n_{ij}$¶
In [27]:
#--- histogram
nn=np.log10(Nij_rnd['R*']*Nij_rnd['T*'])
nmin=nn.min()
nmax=nn.max()
bins=np.linspace(nmin,nmax,int(nmax-nmin)*6)
hist, edge = np.histogram(nn,bins=bins,normed=True)
count, junk = np.histogram(nn,bins=bins)
hist=hist[count!=0]
edge=edge[:-1][count!=0]
count = count[ count !=0 ]
err_hist = hist / count**.5
#--- accumulated histogram
slist=np.array(nn)
slist.sort()
N = len(slist)
d = histogramACCUMLTD( slist.tolist() )
keys=d.keys()
keys.sort()
xx=[];yy=[]
for ikey in keys:
xx.append(d[ikey][0])
yy.append(d[ikey][2])
#------------
#--- plot
#------------
fig= plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.set_xlabel('$log~n_{ij}$')
ax.set_ylabel('$P(log~n_{ij})$')
(nmin,nmax) = (min(slist),max(slist))
ax.set_xlim(floor(nmin),ceil(nmax))
ax.set_xlim(-7,1)
ax.set_ylim(1e-2,1)
#ax.set_xscale('log')
ax.set_yscale('log')
ax.errorbar(edge,hist,yerr=err_hist,color='black',
fmt='-o', markersize=8,markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1)
# ax2 = ax.twinx()
# ax2.plot(xx,yy,
# linestyle='-', drawstyle='steps-post',color='red',
# linewidth=1.0) #--- accumulated
# ax2.set_xlim(-7,1)
# ax2.set_ylim(0,1200)
# ax2.tick_params(axis='y',colors='red')
# ax2.set_ylabel('$N(<n_{ij})$',color='red')
np.save('random',np.c_[edge,hist,err_hist])
#np.save('acc_random',np.c_[xx,yy])
#fig.savefig('%s/pn_shuffled.png'%DIR_OUTPT_figs,bbox_inches='tight')
plt.show()
In [28]:
act = np.load('actual.npy')
rnd = np.load('random.npy')
# acc_act = np.load('acc_actual.npy')
# acc_rnd = np.load('acc_random.npy')
ax.errorbar(edge,hist,yerr=err_hist,color='black',
fmt='-o', markersize=8,markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1)
#------------
#--- plot
#------------
fig= plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.errorbar(act[:,0],act[:,1],yerr=act[:,2],
fmt='-o',color='black',label='Actual',
markersize=8,markerfacecolor='white',markeredgewidth=0.7,
linewidth=.5, barsabove=None,capsize=5,capthick=1,elinewidth=1)
ax.errorbar(rnd[:,0],rnd[:,1],yerr=rnd[:,2],
fmt='-s',color='red',label='Shuffled',
markersize=8,markerfacecolor=None,markeredgewidth=0.7,
linewidth=.5, barsabove=None,capsize=5,capthick=1,elinewidth=1)
(nmin,nmax) = (min(slist),max(slist))
#
ax.set_xlim(-7,1)
ax.set_ylim(1e-2,1)
ax.set_yscale('log')
ax.tick_params(axis='x',top=True)
ax.tick_params(axis='y',right=True,which='both')
#nthh=-2.9
ax.plot([np.log10(n_thresh),np.log10(n_thresh)],[1e-2,0.5])
# ax2 = ax.twinx()
# ax2.plot(acc_act[:,0],acc_act[:,1],
# linestyle='-', drawstyle='steps-post',color='black',
# linewidth=1.0) #--- accumulated
# ax2.plot(acc_rnd[:,0],acc_rnd[:,1],
# linestyle='-', drawstyle='steps-post',color='red',
# linewidth=1.0) #--- accumulated
# ax2.set_xlim(-7,1)
# ax2.set_ylim(0,1200)
# ax2.tick_params(axis='y',colors='red')
# ax2.set_ylabel('$N(<n_{ij})$',color='red')
ax_labels = None #True #False
if ax_labels:
ax.set_xlabel('$log~n_{ij}$')
ax.set_ylabel('$P(log~n_{ij})$')
quantile=5
string = r'HypoDD \newline'
string += r'$m_c=%s,~b=%s,~d_f=%s,~q=%s \newline$'%(mc,bval,Df,int(100*quantile)) #'HypoDD'
string += r'$%s-%s-%s<t<%s-%s-%s$'%(xlo.year,xlo.month, xlo.day, xhi.year,xhi.month, xhi.day) #'HypoDD'
ax.set_title(string,fontsize=12) #'x-small')
ax.legend(fontsize='xx-small')
#DrawFrame(ax, (0.25,0.08),(0.16,0.08),0.1)
DrawFrame(ax, (0.25,0.08),(0.16,0.08),0.1,LOG_Y=True)
fig.savefig('%s/pn_actual_shuffled.png'%DIR_OUTPT_figs,bbox_inches='tight',dpi=2*75)
plt.show()
plot density map¶
In [29]:
#--------------------
#--- plot scattered
#--------------------
fig = plt.figure(figsize=(4,4))
#--------------------
#--- interpolated data
#--------------------
#for nn, index, title in zip([ nij, Nij_rnd ],xrange(2),['actual','randomized']):
for nn, index, title in zip([ nij ],xrange(2),['actual']):
ax = fig.add_subplot(111) #,2,index+1)
#--- interpolate
heatmap=GetInterpolatedData(nn,(TMIN,TMAX),(RMIN,RMAX))
im = ax.pcolormesh(heatmap.T,cmap='coolwarm')
#---
# ax.set_xlabel('T*')
# ax.set_ylabel('R*')
# ax.set_title(title)
##--- reset x-ticks
reset_xticks(4,ax,(TMIN,TMAX))
#--- reset y-ticks
reset_yticks(4,ax,(RMIN,RMAX))
#--- nij must be within limits
nij_copy = nn.copy()
nij_copy=nij_copy[(nij_copy['T*'] >= 10**TMIN ) &
(nij_copy['T*'] < 10**TMAX ) &
(nij_copy['R*'] >= 10**RMIN ) &
(nij_copy['R*'] < 10**RMAX )]
#--- mapping
xx=ax.axis()[0]+(ax.axis()[1]-ax.axis()[0])*(np.log10(nij_copy['T*'])-TMIN)/(TMAX-TMIN)
yy=ax.axis()[2]+(ax.axis()[3]-ax.axis()[2])*(np.log10(nij_copy['R*'])-RMIN)/(RMAX-RMIN)
ax.set_xlim(ax.axis()[0],ax.axis()[1])
ax.set_ylim(ax.axis()[2],ax.axis()[3])
#--- plot
DrawFrame(ax, (0.25,0.08),(0.16,0.08),0.01)
ax.scatter(xx,
yy,
alpha=1,color='black',s=0.01)
#--- draw the separation line
xx=ax.axis()[0]+(ax.axis()[1]-ax.axis()[0])*(np.log10(nij_copy['T*'])-TMIN)/(TMAX-TMIN)
yy=ax.axis()[2]+(ax.axis()[3]-ax.axis()[2])*(np.log10(n_thresh/nij_copy['T*'])-RMIN)/(RMAX-RMIN)
ax.plot( xx, yy,'-.',color='white',linewidth=1)
#fig.colorbar(im,ax=ax)
fig.savefig('%s/densityMaps.png'%DIR_OUTPT_figs,bbox_inches='tight',dpi=150)
ax.set_title('mc=%s,~bval=%s'%(mc,bval),fontsize='xx-small')
#---
mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings
rcc('text', usetex=True)
font = {'size' : 20}
matplotlib.rc('font', **font)
plt.show()
In [30]:
#--- triggering part
nij_trig=nij[nij['T*']*nij['R*']<=n_thresh]
display(nij_trig.info())
trig_ratio = cdf_actual[edges_actual[:-1]<=n_thresh][-1]
print 'trig_ratio=%s'%trig_ratio
In [31]:
#--- plot clusters
def Inside(t,(tmin,tmax)):
if tmin<= t<tmax:
return True
return False
def GetTrigList( nij_trig ):
d_trig = {}
for items in nij_trig.itertuples():
triggerID = items.Parent_id
d_trig.setdefault( triggerID, [] ).append( items.Event_id )
return d_trig
def PlotArrows(df_complete,d_trig, (tmin,tmax), (xmin,xmax), (ymin,ymax) ):
tlist=[]
mlist=[]
for triggerID in d_trig:
t0 = df_complete['date'].iloc[ triggerID ]
if not ( Inside(t0,(tmin,tmax) ) and
Inside(df_complete['longitude'].iloc[triggerID],(xmin,xmax)) and
Inside(df_complete['latitude'].iloc[triggerID],(ymin,ymax)) ):
continue
tlist.append( t0 )
x0 = md.date2num(t0)
y0 = df_complete['magnitude'].iloc[triggerID]
mlist.append( y0 )
for daughter_id in d_trig[triggerID]:
t1 = df_complete['date'].iloc[ daughter_id ]
tlist.append( t1 )
mlist.append( df_complete['magnitude'].iloc[daughter_id] )
if not ( Inside(t1,(tmin,tmax) ) and
Inside(df_complete['longitude'].iloc[daughter_id],(xmin,xmax)) and
Inside(df_complete['latitude'].iloc[daughter_id],(ymin,ymax)) ):
continue
xw = md.date2num(t1) - x0
yw = df_complete['magnitude'].iloc[daughter_id] - y0
plt.annotate("", (t1,df_complete['magnitude'].iloc[daughter_id]), xytext=(t0, y0),
textcoords='data',
arrowprops=dict(arrowstyle="-|>,head_width=.4,head_length=.8",color="b",linewidth="0.3"))
#--- plot circles
df=pd.DataFrame({'date':tlist,'mag':mlist})
plt.scatter(df['date'], df['mag'],
s=4**(df['mag']),
alpha=1,
facecolors='red',color='black')
##plt.savefig('timeSeries.png')
#--------------------------------------
#----- key: event value: aftershock id
#--------------------------------------
d_trig = GetTrigList( nij_trig )
#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi.reset_index(inplace=True,drop=True)
#--------------------
#--- cartesian coords
#--------------------
#swarm_lohi = GetCartesian( swarm_lohi )
#--- setup
fig = plt.figure(figsize=(16,4),frameon=False, dpi=300)
#--- xlimit
xmin = swarm_lohi['longitude'].min()
xmax = swarm_lohi['longitude'].max()
ymin = swarm_lohi['latitude'].min()
ymax = swarm_lohi['latitude'].max()
plt.tick_params(axis='x', labelsize=14)
plt.tick_params(axis='y', labelsize=14)
#--- plot all events
plt.ylim(mc,swarm_lohi['magnitude'].max())
plt.scatter(swarm_lohi['date'],swarm_lohi['magnitude'],
s=2*np.exp(swarm_lohi['magnitude']),
alpha=0.04)
#--- plot arrows
tt=swarm_lohi[swarm_lohi['magnitude']==swarm_lohi['magnitude'].max()]['date'].iloc[0]
PlotArrows( swarm_lohi, d_trig,
( tt+datetime.timedelta(days=-1),tt+datetime.timedelta(days=+1)),
(xmin,xmax),
(ymin,ymax) )
plt.xlim(tt+datetime.timedelta(days=-1),tt+datetime.timedelta(days=+1))
Out[31]:
In [32]:
#--- spatial clusters
# def GetCartesian2nd( dff ):
# df = dff.copy()
# xlo = df['longitude'].min()
# xhi = df['longitude'].max()
# ylo = df['latitude'].min()
# yhi = df['latitude'].max()
# getDistX = lambda x: geopy.distance.vincenty( ( 0.0, xlo ), ( 0.0, x ) ).km
# getDistY = lambda y: geopy.distance.vincenty( ( ylo, 0.0 ), ( y, 0.0 ) ).km
# df[ 'x' ] = df[ 'longitude' ].apply( getDistX )
# df[ 'y' ] = df[ 'latitude' ].apply( getDistY )
# return df
def PlotMap(df_complete,TRIG_ID, (tmin,tmax), (xmin,xmax), (ymin,ymax) ):
#--- plot
counter = 0
exponent=4
for triggerID,colin in zip(TRIG_ID,colors):
t0 = df_complete['date'].iloc[ triggerID ]
# if not ( Inside(t0,(tmin,tmax) ) and
# Inside(df_complete['y(m)'].iloc[triggerID],(xmin,xmax)) and
# Inside(df_complete['z(m)'].iloc[triggerID],(ymin,ymax)) ):
# continue
counter += 1
if counter > 10:
break
#--- plot
fig = plt.figure(figsize=(5,5))#,dpi=150)
ax = fig.add_subplot(111)
# plt.xlim(xmin,xmax)
# plt.ylim(ymin,ymax)
# pdb.set_trace()
print 'm_ms=', df_complete.loc[triggerID]['magnitude']
plt.tick_params(axis='x', labelsize=14)
plt.tick_params(axis='y', labelsize=14)
plt.scatter(df_complete.loc[triggerID]['y(m)'],
df_complete.loc[triggerID]['z(m)'],
s=exponent**df_complete.loc[triggerID]['magnitude'],
facecolors='white',color='white')
d={'x':[],'y':[],'m':[],'t':[]}
#--- plot daughters
for daughter_id in d_trig[triggerID]:
t1 = df_complete['date'].iloc[ daughter_id ]
# if not ( Inside(t1,(tmin,tmax) ) and
# Inside(df_complete['y(m)'].iloc[daughter_id],(xmin,xmax)) and
# Inside(df_complete['z(m)'].iloc[daughter_id],(ymin,ymax)) ):
# continue
#--- remove pairs with r > cs.t
# if rmat[ triggerID, daughter_id ] > cs * tmat[ triggerID, daughter_id ]:
# print 'hello!'
# continue
d['t'].append( t1 )
d['x'].append(df_complete['y(m)'].iloc[daughter_id])
d['y'].append(df_complete['z(m)'].iloc[daughter_id])
d['m'].append(df_complete['magnitude'].iloc[daughter_id])
df=pd.DataFrame(d)
plt.scatter(df['x'],df['y'],
s=exponent**df['m'],
c=df['t'], cmap='jet') #,alpha=0.1)
plt.colorbar()
#--- plot mother
plt.scatter(df_complete['y(m)'].iloc[triggerID],df_complete['z(m)'].iloc[triggerID],
s=exponent**df_complete['magnitude'].iloc[triggerID],
marker='*',
facecolors='yellow',color='black',
alpha=1.0)
#---
colors = itertools.cycle(["r", "b", "g","b"])
#cs = 3 * (24*3600) #--- km per day
#--- sort
new_list = [[len(d_trig[triggerID]),triggerID] for triggerID in d_trig]
new_list.sort(reverse=True)
TRIG_ID = [i[1] for i in new_list]
#--- xlimit
swarm_lohi=GetCartesian(swarm_lohi)
swarm_lohi['y(m)']=swarm_lohi['x(km)']
swarm_lohi['z(m)']=swarm_lohi['z(km)']
xmin = swarm_lohi['y(m)'].min()
xmax = swarm_lohi['y(m)'].max()
ymin = swarm_lohi['z(m)'].min()
ymax = swarm_lohi['z(m)'].max()
#plt.tick_params(axis='x', labelsize=14)
#plt.tick_params(axis='y', labelsize=14)
#--- plot arrows
PlotMap( swarm_lohi, TRIG_ID,
# (pd.to_datetime('2010-04-04'),pd.to_datetime('2010-06-26')),
# (pd.to_datetime('2010-04-08'),pd.to_datetime('2010-04-14')),
( swarm['date'].min(), swarm['date'].max()),
(xmin,xmax),
(ymin,ymax) )
In [33]:
#--- cluster life time
def GetPairsWithSpecifiedParentMag((m0,m1),catalog,df):
df_parents = df.groupby(by='Parent_id').groups #--- parent events
ds = catalog['magnitude'].loc[df_parents.keys()] #--- mag. of parent events
ds_m=ds[(m0 <= ds) & (ds < m1)] #--- parent events with m0<m<m1
df_m=df[pd.DataFrame([df['Parent_id'] ==k for k in ds_m.index]).any()] #--- data frame
return df_m
#def AddDist( df_trig, df_complete ):
# x=df_complete.loc[ df_trig['Event_id'] ]['r(m)']
# y=df_complete.loc[ df_trig['Parent_id'] ]['r(m)']
# df_trig['R'] = pd.Series(numpy.linalg.norm(x-y)))
# assert len ( df_trig[ df_trig['R'] == 0.0 ] ) == 0, '%s'%display( df_trig[ df_trig['R'] == 0.0 ] )
def AddTime( df_trig, df_complete ):
prefact = 1.0 / ( 24.0 * 60 * 60 ) #--- daily
x=df_complete.loc[ df_trig['Event_id'] ]
y=df_complete.loc[ df_trig['Parent_id'] ]
x.reset_index(inplace=True)
y.reset_index(inplace=True)
df_trig['T'] = x['date']-y['date']
df_trig['T'] = df_trig['T'].apply(lambda x: x.total_seconds() * prefact)
def RemovePair( df, cs ):
return df [ df[ 'R' ] <= df[ 'T' ] * cs ]
#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi.reset_index(inplace=True,drop=True)
#--- choose m0<m<m1
n = 7
#m = list(np.linspace(mc,swarm['magnitude'].max(),n))
m = list(np.linspace(mc,swarm_lohi['magnitude'].max()+1e-6,n))
#m+=[swarm['magnitude'].max()]
X={}
rho_dict={}
m_range={}
Err={}
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
for i in xrange( len(m) - 1 ):
m0 = m[ i ]
m1 = m[ i + 1 ]
# print m0, m1
# nij_trig.groupby(by='Parent_id')
df_trig = GetPairsWithSpecifiedParentMag((m0,m1),swarm_lohi,nij_trig) #--- get parent with m0<m<m1
df_trig.reset_index( inplace = True, drop = True )
#--- add distance & time column
# df_trig['R'] = rmat[df_trig['Parent_id'],df_trig['Event_id']]
# df_trig['T'] = tmat[df_trig['Parent_id'],df_trig['Event_id']]
# assert len ( df_trig[ df_trig['R'] == 0.0 ] ) == 0, '%s'%display( df_trig[ df_trig['R'] == 0.0 ] )
# assert len ( df_trig[ df_trig['T'] == 0.0 ] ) == 0, '%s'%display( df_trig[ df_trig['T'] == 0.0 ] )
#--- remove pairs with r > cs.t
# cs = 3 * (24*3600) #--- km per day
# df_trig = RemovePair( df_trig, cs )
#--- life time
df_duration = df_trig.groupby(by='Parent_id')['T'].max()
#--- rho plots
n_decades = int( math.ceil( log(df_duration.max()/df_duration.min(),10) ) )
nbin_per_decade = 3
nbins = nbin_per_decade * n_decades
rho, xedges=np.histogram(df_duration,
density=True,
bins=np.logspace(log(df_duration.min(),10),
log(df_duration.max(),10),nbins))
hist, xedges=np.histogram(df_duration,
bins=np.logspace(log(df_duration.min(),10),
log(df_duration.max(),10),nbins))
err = rho/np.sqrt(hist)
#--- rho
# plt.xlim(1e-3,1e2)
# plt.ylim(1e-3,1e2)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Duration(Day)')
plt.ylabel('P(Duration)')
plt.errorbar(xedges[:-1],rho, yerr=err,label='[%2.1f,%2.1f]'%(m0,m1))
fontP = FontProperties()
fontP.set_size('small')
plt.legend(title='m',prop=fontP,loc='upper right',bbox_to_anchor=(1.5, 1))
x=xedges[:-1]
X[i] = x[:]
rho_dict[i]=rho[:]
m_range[i]=np.array([m0,m1])
Err[i]=err[:]
#plt.savefig('%s/duration.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight')
In [34]:
#--- histograms: magnitude of mother events
#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi.reset_index(inplace=True,drop=True)
mlist = swarm_lohi.loc[ nij_trig.groupby(by='Parent_id').groups.keys() ]['magnitude']
#plt.hist(mlist)
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
mlist0= swarm_lohi['magnitude'].loc[df_zeros.index]
#--- accumulated histogram
N = len(mlist)+len(mlist0)
slist=mlist.tolist()+mlist0.tolist() #np.array(swarm_copy['magnitude'])
slist.sort()
d = histogramACCUMLTD( slist )
keys=d.keys()
keys.sort()
#--- plot
fig, ax = plt.subplots(figsize=(4,4))
ax.set_yscale('log')
ax.set_xlabel('$m$')
ax.set_ylabel('$N(m)$')
ax.axis([1,6,1,10000])
ax.set_xticks(np.linspace(1,6,6))
xmin=min(slist)
xmax=max(slist)
dx=0.1
junk = ax.hist( slist,
bins=int((xmax-xmin)/dx),
label='histogram',color='silver') #--- histogram
xx=[];yy=[]
for ikey in keys:
xx.append(d[ikey][0])
yy.append(N-d[ikey][2])
ax.plot(xx,yy,
linestyle='-', drawstyle='steps-post',color='black',
linewidth=1.0) #--- accumulated
b=0.95
ax.plot([0,6],[1e6, 1e6* 10**(-b*(6-0))],'r-.',linewidth=1)
#DrawFrame(ax, (0.2,0.06),(0.14,0.06),0.1,LOG_Y=True)
#plt.savefig('%s/gr.png'%DIR_OUTPT_figs,dpi=300,bbox_inches='tight',transparent=True)
Out[34]:
In [35]:
#--- productivity relation
import traceback
mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings
rcc('text', usetex=True)
font = {'size' : 26}
matplotlib.rc('font', **font)
import numpy as np
import math
def GetPairsWithSpecifiedParentMag((m0,m1),catalog,df):
df_parents = df.groupby(by='Parent_id').groups #--- parent events
ds = catalog['magnitude'].loc[df_parents.keys()] #--- mag. of parent events
ds_m=ds[(m0 <= ds) & (ds < m1)] #--- parent events with m0<m<m1
df_m=df[pd.DataFrame([df['Parent_id'] ==k for k in ds_m.index]).any()] #--- data frame
return df_m, np.mean( ds_m ), np.std( ds_m )/(len(ds_m))**0.5 #--- exponential: std=mean
#def AddDist( df_trig, df_complete ):
# x=df_complete.loc[ df_trig['Event_id'] ]['r(m)']
# y=df_complete.loc[ df_trig['Parent_id'] ]['r(m)']
# df_trig['R'] = pd.Series(numpy.linalg.norm(x-y)))
# assert len ( df_trig[ df_trig['R'] == 0.0 ] ) == 0, '%s'%display( df_trig[ df_trig['R'] == 0.0 ] )
def AddTime( df_trig, df_complete ):
prefact = 1.0 / ( 24.0 * 60 * 60 ) #--- daily
x=df_complete.loc[ df_trig['Event_id'] ]
y=df_complete.loc[ df_trig['Parent_id'] ]
x.reset_index(inplace=True)
y.reset_index(inplace=True)
df_trig['T'] = x['date']-y['date']
df_trig['T'] = df_trig['T'].apply(lambda x: x.total_seconds() * prefact)
def RemovePair( dff, cs ):
df = dff.copy()
df['R'] = rmat[df['Parent_id'],df['Event_id']]
df['T'] = tmat[df['Parent_id'],df['Event_id']]
assert len ( df[ df['R'] == 0.0 ] ) == 0, '%s'%display( df[ df['R'] == 0.0 ] )
assert len ( df[ df['T'] == 0.0 ] ) == 0, '%s'%display( df[ df['T'] == 0.0 ] )
return df [ df[ 'R' ] <= df[ 'T' ] * cs ]
def AdjustBins(data,bins):
hist, xedges=np.histogram(data, #--- histogram
bins=bins)
tmp = [xedges[0]]
for i,j in zip(xedges[1:],hist): #--- expand bins
if j < 2:
continue
tmp.append(i)
# tmp.append(xedges[-1])
if hist[-1] < 2:
tmp.pop(-2)
return tmp
def GetParentsMagnitudes( swarm_lohi, df_trig ):
mlist = swarm_lohi.loc[ df_trig.groupby(by='Parent_id').groups.keys() ]['magnitude']
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
mlist0= swarm_lohi['magnitude'].loc[df_zeros.index]
return mlist.tolist()+mlist0.tolist()
#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi.reset_index(inplace=True,drop=True)
#--- remove pairs with r > cs.t
#cs = 3 * (24*3600) #--- km per day
df_trig = nij_trig.copy()
#--- parents' magnitudes
m_list = GetParentsMagnitudes( swarm_lohi, df_trig )
#--- number of daughters for each mother
df_n = df_trig.groupby(by='Parent_id').count()['Event_id']
#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
#--- scattered plot
fig = plt.figure(figsize=(4,4),frameon=True)
ax = fig.add_subplot(111)
m_min = min(m_list)-1.0e-6
m_max = max(m_list)+1.0e-6 #--- exclude the largest event
ax.set_xlim(np.floor(m_min),np.ceil(m_max))
ax.set_ylim(0.1,1000)
ax.plot( swarm_lohi['magnitude'].loc[df_n.index],
df_n,
'x',color='black',alpha=0.5,zorder=1,label='')
ax.legend(fontsize=0)
ax2=ax.twinx()
ax2.set_yscale('linear')
ax2.set_xlim(np.floor(m_min),np.ceil(m_max))
ax2.set_ylim(-.02,1)
ax2.plot( swarm_lohi['magnitude'].loc[df_zeros.index],
df_zeros,
'x',color='black',alpha=0.25)
ax2.tick_params(labelright=False,labelleft=False,left=False, right=False)
ax.set_yscale('log')
ax.set_ylim(.11,1000)
#plt.xlabel('$m$')
#plt.ylabel(r'$N_{as}$')
ax.tick_params(axis='y',left=True, right=True,which='both')
ax.tick_params(axis='x',bottom=True, top=True,which='both')
#--- choose m0<m<m1
dm = 0.1
n = int( ( m_max - m_min ) / dm )
m = list(np.linspace(m_min,m_max,n))
m=AdjustBins(m_list,m)
Mlist=m
m_mean=np.array(np.zeros(len(m) - 1))
m_err=np.array(np.zeros(len(m) - 1))
N_as = np.array(np.zeros(len(m) - 1))
err_N = np.array(np.zeros(len(m) - 1))
for i in xrange( len(m) - 1 ):
try:
m0 = m[ i ]
m1 = m[ i + 1 ]
# print m0, m1
df, m_mean[ i ], m_err[i] = GetPairsWithSpecifiedParentMag((m0,m1),
swarm_lohi,df_trig) #--- get parent with m0<m<m1
df.reset_index( inplace = True, drop = True )
#--- number of daughters for each mother
df_n = df.groupby(by='Parent_id').count()['Event_id']
#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
subset = swarm_lohi.loc[list_of_zero_afs]
list_of_zero_afs = subset[(subset['magnitude'] >= m0) & (subset['magnitude'] <= m1)].index
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
#--- mean
N_as[ i ] = 1.0*(df_n.sum()+ df_zeros.sum())/(len(df_n)+len(df_zeros))
err_N[ i ] = (N_as[ i ]/(len(df_n)+len(df_zeros)))**0.5 #--- poisson
# print N_as[ i ]
# err_N[ i ]=np.std(df_n.tolist()+df_zeros.tolist())/(len(df_n)+len(df_zeros))**0.5
except:
# traceback.print_exc()
continue
#ax.plot(m, np.append(N_as,N_as[-1]), linestyle='-', drawstyle='steps-post',color='red')
ax.errorbar(m_mean,N_as,
yerr=err_N,xerr=m_err,
markerfacecolor='white',
fmt='o', markersize=8,color='red',markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1)
#--- add major xticks
xmin=np.ceil(ax.axis()[0])
xmax=np.floor(ax.axis()[1])
nbin = xmax - xmin
ax.set_xticks(np.linspace(xmin,xmax,int(nbin)+1))
DrawFrame(ax, (0.24,0.08),(0.14,0.06),0.01,LOG_Y=True)
#from scipy.optimize import curve_fit
def func(x,alpha,c):
return c*10**(alpha*x)
#popt, pconv = curve_fit(func,m_mean,N_as)
popt = np.polyfit(m_mean,np.log10(N_as), 1)
popt[1] = 10**(popt[1])
alpha_p = popt[0]
mm = swarm_lohi['magnitude'].loc[df_zeros.index].tolist() + \
swarm_lohi['magnitude'].loc[df_n.index].tolist() #swarm_lohi[swarm_lohi['magnitude']>=mc]['magnitude']
mm = np.array(mm)
mm=list(set(mm))
mm.sort()
mm=np.array(mm)
ax.plot(mm,func(mm,*popt),'-.g',zorder=3,label='$N\propto 10^{%3.2f m}$'%popt[0])
ax.legend(fontsize='xx-small',frameon=False)
#ax.set_title('q=%d%%,df=%s'%(quantile*100,Df),fontsize='xx-small')
#print 'alpha=%s'%popt[0]
#print pconv**0.5
plt.savefig('%s/prod.png'%DIR_OUTPT_figs,dpi=75*.75,bbox_inches='tight')
np.save('%s/alpha'%DIR_OUTPT_figs,alpha_p)
mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings
rcc('text', usetex=True)
font = {'size' : 20 }
matplotlib.rc('font', **font)
plt.show()
based on quantiles¶
In [36]:
#--- productivity relation
mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings
rcc('text', usetex=True)
font = {'size' : 26}
matplotlib.rc('font', **font)
import numpy as np
import math
from astropy import stats
from scipy.stats import chi2
def AdjustBins(data,bins):
hist, xedges=np.histogram(data, #--- histogram
bins=bins)
tmp = [xedges[0]]
for i,j in zip(xedges[1:],hist): #--- expand bins
if j < 2:
continue
tmp.append(i)
tmp.append(xedges[-1])
# if hist[-1] < 2:
# tmp.pop(-2)
return tmp
def GetEvenQuantiles(mlist, nn):
N=len(mlist)
hist, edges = np.histogram(mlist,bins=1024*10,density=None)
assert np.sum(hist) == N
# plt.plot(edges[:-1],np.cumsum(hist))
# plt.show()
Quantiles = np.concatenate([(1.0*np.arange(nn)/nn)[1:],np.array([1])])
# print Quantiles
m_q=[np.min(mlist)-1e-6]
for q in Quantiles:
m_q+=[edges[:-1][np.cumsum(hist)<=q*N][-1]]
# print len(mlist[mlist<m_q[1]]),N/2
return m_q
def GetPval( m_mean, N_as, err_N, func, popt ):
x = m_mean
y = N_as
ye = err_N
#--- goodness of fit
n=len(x)
m=2 #--- number of coefficients
xsq = np.sum(((y-func(x,*popt))/ye)**2)
pval=1-chi2.cdf(xsq,n-m)
return pval
def PointWiseVariance( m_mean, N_as):
n = len( m_mean )
p = 1
y = np.log10(N_as)
x = m_mean
#--- estimate beta
X = np.c_[x,np.ones( n )]
beta = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T,X)),X.T),y)
#--- estimate sigma
def func( x, beta0, beta1 ):
return beta0 * x + beta1
y_hat = func( x, *beta ) #--- estimated output
sigma_sq = np.sum( ( y - y_hat ) ** 2 ) / ( n - p - 1 ) #--- variance
var_beta = np.linalg.inv(np.matmul(X.T,X)) * sigma_sq #--- uncertainties in beta
var_y = np.matmul( np.matmul( X, var_beta ), X.T ) #--- pointwise
y_up = y_hat + 2 * np.diag(var_y) ** 0.5 #--- std-error bands
y_lo = y_hat - 2 * np.diag(var_y) ** 0.5
return y_lo, y_up
def ProdRelation( nij_trig, swarm_lohi, **kwargs ):
# print kwargs
#--- remove pairs with r > cs.t
df_trig = nij_trig.copy()
#--- parents' magnitudes
m_list = GetParentsMagnitudes( swarm_lohi, df_trig )
#--- number of daughters for each mother
df_nn = df_trig.groupby(by='Parent_id').count()['Event_id']
#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
# print swarm_lohi['magnitude'].loc[df_zeros.index]
# print df_zeros
#--- parents' magnitudes including those with zero child
mm = swarm_lohi['magnitude'].loc[df_zeros.index].tolist() + \
swarm_lohi['magnitude'].loc[df_nn.index].tolist() #swarm_lohi[swarm_lohi['magnitude']>=mc]['magnitude']
mm = np.array(mm)
# mm=list(set(mm))
mm.sort()
mlist_with_zero=np.array(mm)
#--- choose m0<m<m1
dm = 0.1
m_min = min(mlist_with_zero)-1.0e-6 #0.5*dm
m_max = max(mlist_with_zero)+1.0e-6 #0.5*dm
m_lo = kwargs['mlo'] if 'mlo' in kwargs else m_min
m_hi = kwargs['mhi'] if 'mhi' in kwargs else m_max
n = int( ( m_max - m_min ) / dm )
m = list(np.linspace(m_min,m_max,n))
#
Mlist = m=np.histogram_bin_edges(mlist_with_zero, bins='auto')
if 'NON_UNIFRM' in kwargs and kwargs['NON_UNIFRM']:
Mlist = m = stats.bayesian_blocks(mlist_with_zero)
# dxx, m = stats.knuth_bin_width(mlist_with_zero, return_bins=True)
#
if 'QUANT' in kwargs and kwargs['QUANT']:
Mlist = m = GetEvenQuantiles(mlist_with_zero, n)
#
if 'ADJ_BIN' in kwargs and kwargs[ 'ADJ_BIN' ]:
m=AdjustBins(mlist_with_zero,m)#[:-1]
if 'MBINS' in kwargs and kwargs[ 'MBINS' ]:
m = kwargs[ 'MBINS' ]
Mlist=m
# print len(m)
#--- choose mag. range based on evenly dist. q
# print m, np.max(mlist_with_zero)
m_mean=np.array(np.zeros(len(m) - 1))
m_err=np.array(np.zeros(len(m) - 1))
N_as = np.array(np.zeros(len(m) - 1))
err_N = np.array(np.zeros(len(m) - 1))
for i in xrange( len(m) - 1 ):
try:
m0 = m[ i ] - 1.0e-6
m1 = m[ i + 1 ] + 1.0e-6
df, m_mean[ i ], m_err[i] = GetPairsWithSpecifiedParentMag((m0,m1),
swarm_lohi,df_trig) #--- get parent with m0<m<m1
df.reset_index( inplace = True, drop = True )
#--- number of daughters for each mother
df_n = df.groupby(by='Parent_id').count()['Event_id']
#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
subset = swarm_lohi.loc[list_of_zero_afs]
list_of_zero_afs = subset[(subset['magnitude'] >= m0) & (subset['magnitude'] <= m1)].index
df_zeross = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
N_as[ i ] = 1.0*(df_n.sum()+ df_zeross.sum())/(len(df_n)+len(df_zeross))
err_N[ i ] = (N_as[ i ]/(len(df_n)+len(df_zeross)))**0.5 #--- poisson
# err_N[ i ]=np.std(df_n.tolist()+df_zeros.tolist())/(len(df_n)+len(df_zeros))**0.5
# print len(df_n)+len(df_zeross)
# display(df_n.head())
except:
# print traceback.print_exc()
continue
#--- scattered plot
fig = plt.figure(figsize=(4,4),frameon=True)
ax = fig.add_subplot(111)
title = kwargs['title'] if 'title' in kwargs else ''
ax.set_title(title,fontsize='xx-small')
ax.set_xlim(np.floor(m_min),np.ceil(swarm_lohi['magnitude'].max()))
ax.set_ylim(0.11,10**np.ceil(np.log10(df_nn.max())))
ax.plot( swarm_lohi['magnitude'].loc[df_nn.index],
df_nn,
'x',color='black',alpha=0.5,zorder=2,label='')
ax.plot(m, np.append(N_as,N_as[-1]), linestyle='-', drawstyle='steps-post',color='red',zorder=1)
ax.legend(fontsize=0)
#
ax2=ax.twinx()
ax2.set_yscale('linear')
ax2.set_xlim(np.floor(m_min),np.ceil(swarm_lohi['magnitude'].max()))
ax2.set_ylim(-.02,1)
ax2.plot( swarm_lohi['magnitude'].loc[df_zeros.index],
df_zeros,
'x',color='black',alpha=0.25)
ax2.tick_params(labelright=False,labelleft=False,left=False, right=False)
#
ax.set_yscale('log')
#plt.xlabel('$m$')
#plt.ylabel(r'$N_{as}$')
ax.tick_params(axis='y',left=True, right=True,which='both')
ax.tick_params(axis='x',bottom=True, top=True,which='both')
#--- plot binned data
ax.errorbar(m_mean,N_as,
yerr=err_N,xerr=m_err,
markerfacecolor='white',
fmt='o', markersize=8,color='red',markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1,
zorder=3)
ax.set_xlim(0,6)
ax.set_ylim(.11,1e3)
#--- add major xticks
xmin=np.ceil(ax.axis()[0])
xmax=np.floor(ax.axis()[1])
nbin = xmax - xmin
ax.set_xticks(np.linspace(xmin,xmax,int(nbin)+1))
DrawFrame(ax, (0.24,0.08),(0.14,0.06),0.01,LOG_Y=True)
#--- least-squares fits
def func(x,alpha,c):
return c*10**(alpha*x)
#
#--- filter
m_mean = m_mean[N_as != 0.0]
err_N = err_N[N_as != 0.0]
N_as = N_as[N_as != 0.0]
#
indices = np.all([m_mean >= m_lo,m_mean < m_hi],axis=0)
m_mean = m_mean[indices]
err_N = err_N[indices]
N_as = N_as[indices]
#
err_log = err_N / N_as / np.log(10.0)
popt, cov = np.polyfit(m_mean,np.log10(N_as), 1, cov=True, w=1.0 / err_log) #'unscaled')
popt[1] = 10**(popt[1])
alpha_p = popt[0]
#--- pvalue
pval_lin = GetPval( m_mean, N_as, err_N, func, popt )
#--- plot the fit
# ax.plot(m_mean,func(m_mean,*popt),'-.g',zorder=3,label='linear$N\propto 10^{%3.2f m},p=%3.2e$'%(popt[0],pval_lin))
#print 'alpha=%s'%popt[0]
#print pconv**0.5
# #--- point-wise variance
# # y_lo, y_up = PointWiseVariance( m_mean, N_as)
# # #--- plot
# # ax.fill_between(m_mean,
# # 10**y_up,
# # 10**y_lo,
# # color='yellow',zorder=2)
#--- non-linear fits
from scipy import optimize
def func( x, alpha, K ):
return K*10**(alpha*x)
xdata = m_mean
ydata = N_as
popt, pcov = optimize.curve_fit(func, xdata, ydata, sigma=err_N)
print 'covariance',pcov
#--- pvalue
pval_nonlin = GetPval( m_mean, N_as, err_N, func, popt )
#--- plot the fit
ax.plot(m_mean,func(m_mean,*popt),'-.b',zorder=4,label='$N\propto 10^{%3.2f m}, p=%3.2e$'%(popt[0],pval_nonlin))
ax.legend(fontsize=13,frameon=False)
#--- save
dpi = kwargs['dpi'] if 'dpi' in kwargs else 150
plt.savefig('%s/prod_%s.png'%(DIR_OUTPT_figs,title),dpi=dpi,bbox_inches='tight')
mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings
rcc('text', usetex=True)
font = {'size' : 20}
matplotlib.rc('font', **font)
plt.show()
return pval_lin, pval_nonlin, m
#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi.reset_index(inplace=True,drop=True)
pval_lin, pval_nonlin, mbins = ProdRelation( nij_trig, swarm_lohi,
# mlo=mc-0.05, mhi=3.6+0.05,
# MBINS = [0.9296289062500129, 1.0296289062500126, 1.1296289062500122, 1.2296289062500119, 1.3296289062500115, 1.4296289062500112, 1.5296289062500108, 1.6296289062500104, 1.72962890625001, 1.8296289062500097, 1.9296289062500094, 2.029628906250009, 2.1296289062500087, 2.2296289062500083, 2.329628906250008, 2.4296289062500076, 2.5296289062500072, 2.629628906250007, 2.7296289062500065, 2.829628906250006, 2.929628906250006, 3.0296289062500055, 3.129628906250005, 3.2296289062500048, 3.3296289062500044, 3.429628906250004, 3.5296289062500037, 3.6296289062500033, 3.729628906250003, 3.8296289062500026, 3.9296289062500023, 4.029628906250002, 4.129628906250002, 4.229628906250001, 4.329628906250001, 4.4296289062500005, 4.52962890625, 5.72],
NON_UNIFRM = False, QUANT = False, ADJ_BIN = False,
title='Yuha', dpi = 67)
#print mlist
#np.save('%s/alpha'%DIR_OUTPT_figs,alpha_p)
In [37]:
#--- productivity relation
mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings
rcc('text', usetex=True)
font = {'size' : 26}
matplotlib.rc('font', **font)
import numpy as np
import math
from astropy import stats
from scipy.stats import chi2
def AdjustBins(data,bins):
hist, xedges=np.histogram(data, #--- histogram
bins=bins)
tmp = [xedges[0]]
for i,j in zip(xedges[1:],hist): #--- expand bins
if j < 2:
continue
tmp.append(i)
tmp.append(xedges[-1])
# if hist[-1] < 2:
# tmp.pop(-2)
return tmp
def GetEvenQuantiles(mlist, nn):
N=len(mlist)
hist, edges = np.histogram(mlist,bins=1024*10,density=None)
assert np.sum(hist) == N
# plt.plot(edges[:-1],np.cumsum(hist))
# plt.show()
Quantiles = np.concatenate([(1.0*np.arange(nn)/nn)[1:],np.array([1])])
# print Quantiles
m_q=[np.min(mlist)-1e-6]
for q in Quantiles:
m_q+=[edges[:-1][np.cumsum(hist)<=q*N][-1]]
# print len(mlist[mlist<m_q[1]]),N/2
return m_q
def GetPval( m_mean, N_as, err_N, func, popt ):
x = m_mean
y = N_as
ye = err_N
#--- goodness of fit
n=len(x)
m=2 #--- number of coefficients
xsq = np.sum(((y-func(x,*popt))/ye)**2)
pval=1-chi2.cdf(xsq,n-m)
return pval
def PointWiseVariance( m_mean, N_as):
n = len( m_mean )
p = 1
y = np.log10(N_as)
x = m_mean
#--- estimate beta
X = np.c_[x,np.ones( n )]
beta = np.matmul(np.matmul(np.linalg.inv(np.matmul(X.T,X)),X.T),y)
#--- estimate sigma
def func( x, beta0, beta1 ):
return beta0 * x + beta1
y_hat = func( x, *beta ) #--- estimated output
sigma_sq = np.sum( ( y - y_hat ) ** 2 ) / ( n - p - 1 ) #--- variance
var_beta = np.linalg.inv(np.matmul(X.T,X)) * sigma_sq #--- uncertainties in beta
var_y = np.matmul( np.matmul( X, var_beta ), X.T ) #--- pointwise
y_up = y_hat + 2 * np.diag(var_y) ** 0.5 #--- std-error bands
y_lo = y_hat - 2 * np.diag(var_y) ** 0.5
return y_lo, y_up
def ProdRelation( nij_trig, swarm_lohi, **kwargs ):
# print kwargs
#--- remove pairs with r > cs.t
df_trig = nij_trig.copy()
#--- parents' magnitudes
m_list = GetParentsMagnitudes( swarm_lohi, df_trig )
#--- number of daughters for each mother
df_nn = df_trig.groupby(by='Parent_id').count()['Event_id']
#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
# print swarm_lohi['magnitude'].loc[df_zeros.index]
# print df_zeros
#--- parents' magnitudes including those with zero child
mm = swarm_lohi['magnitude'].loc[df_zeros.index].tolist() + \
swarm_lohi['magnitude'].loc[df_nn.index].tolist() #swarm_lohi[swarm_lohi['magnitude']>=mc]['magnitude']
mm = np.array(mm)
# mm=list(set(mm))
mm.sort()
mlist_with_zero=np.array(mm)
#--- choose m0<m<m1
dm = 0.1
m_min = min(mlist_with_zero)-1.0e-6 #0.5*dm
m_max = max(mlist_with_zero)+1.0e-6 #0.5*dm
m_lo = kwargs['mlo'] if 'mlo' in kwargs else m_min
m_hi = kwargs['mhi'] if 'mhi' in kwargs else m_max
n = int( ( m_max - m_min ) / dm )
m = list(np.linspace(m_min,m_max,n))
#
Mlist = m=np.histogram_bin_edges(mlist_with_zero, bins='auto')
if 'NON_UNIFRM' in kwargs and kwargs['NON_UNIFRM']:
Mlist = m = stats.bayesian_blocks(mlist_with_zero)
# dxx, m = stats.knuth_bin_width(mlist_with_zero, return_bins=True)
#
if 'QUANT' in kwargs and kwargs['QUANT']:
Mlist = m = GetEvenQuantiles(mlist_with_zero, n)
#
if 'ADJ_BIN' in kwargs and kwargs[ 'ADJ_BIN' ]:
m=AdjustBins(mlist_with_zero,m)#[:-1]
if 'MBINS' in kwargs and kwargs[ 'MBINS' ]:
m = kwargs[ 'MBINS' ]
Mlist=m
# print len(m)
#--- choose mag. range based on evenly dist. q
# print m, np.max(mlist_with_zero)
m_mean=np.array(np.zeros(len(m) - 1))
m_err=np.array(np.zeros(len(m) - 1))
N_as = np.array(np.zeros(len(m) - 1))
err_N = np.array(np.zeros(len(m) - 1))
for i in xrange( len(m) - 1 ):
try:
m0 = m[ i ] - 1.0e-6
m1 = m[ i + 1 ] + 1.0e-6
df, m_mean[ i ], m_err[i] = GetPairsWithSpecifiedParentMag((m0,m1),
swarm_lohi,df_trig) #--- get parent with m0<m<m1
df.reset_index( inplace = True, drop = True )
#--- number of daughters for each mother
df_n = df.groupby(by='Parent_id').count()['Event_id']
#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
subset = swarm_lohi.loc[list_of_zero_afs]
list_of_zero_afs = subset[(subset['magnitude'] >= m0) & (subset['magnitude'] <= m1)].index
df_zeross = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
N_as[ i ] = 1.0*(df_n.sum()+ df_zeross.sum())/(len(df_n)+len(df_zeross))
err_N[ i ] = (N_as[ i ]/(len(df_n)+len(df_zeross)))**0.5 #--- poisson
# err_N[ i ]=np.std(df_n.tolist()+df_zeros.tolist())/(len(df_n)+len(df_zeros))**0.5
# print len(df_n)+len(df_zeross)
# display(df_n.head())
except:
# print traceback.print_exc()
continue
LABELS = kwargs['LABELS'] if 'LABELS' in kwargs else ''
#--- scattered plot
fig = plt.figure(figsize=(4,4),frameon=True)
ax = fig.add_subplot(111)
title = kwargs['title'] if 'title' in kwargs else ''
if LABELS:
ax.set_title(title,fontsize='xx-small')
ax.set_xlim(np.floor(m_min),np.ceil(swarm_lohi['magnitude'].max()))
ax.set_ylim(0.11,10**np.ceil(np.log10(df_nn.max())))
ax.plot( swarm_lohi['magnitude'].loc[df_nn.index],
df_nn,
'x',color='black',alpha=0.5,zorder=2,label='')
ax.plot(m, np.append(N_as,N_as[-1]), linestyle='-', drawstyle='steps-post',color='red',zorder=1)
ax.legend(fontsize=0)
#
ax2=ax.twinx()
ax2.set_yscale('linear')
ax2.set_xlim(np.floor(m_min),np.ceil(swarm_lohi['magnitude'].max()))
ax2.set_ylim(-.02,1)
ax2.plot( swarm_lohi['magnitude'].loc[df_zeros.index],
df_zeros,
'x',color='black',alpha=0.25)
ax2.tick_params(labelright=False,labelleft=False,left=False, right=False)
#
ax.set_yscale('log')
#plt.xlabel('$m$')
#plt.ylabel(r'$N_{as}$')
ax.tick_params(axis='y',left=True, right=True,which='both')
ax.tick_params(axis='x',bottom=True, top=True,which='both')
#--- plot binned data
ax.errorbar(m_mean,N_as,
yerr=err_N,xerr=m_err,
markerfacecolor='white',
fmt='o', markersize=8,color='red',markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1,
zorder=3)
ax.set_xlim(0,6)
ax.set_ylim(.11,1e3)
#--- add major xticks
xmin=np.ceil(ax.axis()[0])
xmax=np.floor(ax.axis()[1])
nbin = xmax - xmin
ax.set_xticks(np.linspace(xmin,xmax,int(nbin)+1))
DrawFrame(ax, (0.24,0.08),(0.14,0.06),0.1,LOG_Y=True)
#--- least-squares fits
def func(x,alpha,c):
return c*10**(alpha*x)
#
#--- filter
m_mean = m_mean[N_as != 0.0]
err_N = err_N[N_as != 0.0]
N_as = N_as[N_as != 0.0]
#
indices = np.all([m_mean >= m_lo,m_mean < m_hi],axis=0)
m_mean = m_mean[indices]
err_N = err_N[indices]
N_as = N_as[indices]
#
err_log = err_N / N_as / np.log(10.0)
popt, cov = np.polyfit(m_mean,np.log10(N_as), 1, cov=True, w=1.0 / err_log) #'unscaled')
popt[1] = 10**(popt[1])
alpha_p = popt[0]
#--- pvalue
pval_lin = GetPval( m_mean, N_as, err_N, func, popt )
#--- plot the fit
# ax.plot(m_mean,func(m_mean,*popt),'-.g',zorder=3,label='linear$N\propto 10^{%3.2f m},p=%3.2e$'%(popt[0],pval_lin))
#print 'alpha=%s'%popt[0]
#print pconv**0.5
# #--- point-wise variance
# # y_lo, y_up = PointWiseVariance( m_mean, N_as)
# # #--- plot
# # ax.fill_between(m_mean,
# # 10**y_up,
# # 10**y_lo,
# # color='yellow',zorder=2)
#--- non-linear fits
from scipy import optimize
def func( x, alpha, K ):
return K*10**(alpha*x)
xdata = m_mean
ydata = N_as
popt, pcov = optimize.curve_fit(func, xdata, ydata, sigma=err_N)
print 'covariance',pcov
#--- pvalue
pval_nonlin = GetPval( m_mean, N_as, err_N, func, popt )
#--- plot the fit
EXTEND = kwargs['EXTEND'] if 'EXTEND' in kwargs else ''
mrr = np.array( [np.min(mlist_with_zero), np.max(mlist_with_zero)] )
lw = 2.0
string = r'$\newline m\in(%2.1f,%2.1f)$'%(m_lo,m_hi)\
+r'$\newline\newline \alpha=%3.2f$'%popt[0]+' '\
+r'$\newline\newline p$'+r'-val'+'$=%4.3f$'%(pval_nonlin)
if EXTEND:
ax.plot(mrr,func(mrr,*popt),
'-.g',zorder=4, linewidth=lw,
label=string)
else:
ax.plot(m_mean,func(m_mean,*popt),
'-.g',zorder=4, linewidth=lw,
label=string)
if LABELS:
ax.legend(fontsize=19,frameon=False,bbox_to_anchor=(0.57, 0.94),handlelength=0)
#--- save
dpi = kwargs['dpi'] if 'dpi' in kwargs else 150
if 'save' in kwargs:
plt.savefig('%s.png'%kwargs['save'],dpi=dpi,bbox_inches='tight')
mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings
rcc('text', usetex=True)
font = {'size' : 20}
matplotlib.rc('font', **font)
plt.show()
return pval_lin, pval_nonlin, m, popt[ 0 ]
#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi.reset_index(inplace=True,drop=True)
path_save = '.'
pval_lin, pval_nonlin, mbins, alpha_p = ProdRelation( nij_trig, swarm_lohi,
mlo=mc-0.01, mhi=4.4, #3.6+0.05,
# MBINS = [0.9296289062500129, 1.0296289062500126, 1.1296289062500122, 1.2296289062500119, 1.3296289062500115, 1.4296289062500112, 1.5296289062500108, 1.6296289062500104, 1.72962890625001, 1.8296289062500097, 1.9296289062500094, 2.029628906250009, 2.1296289062500087, 2.2296289062500083, 2.329628906250008, 2.4296289062500076, 2.5296289062500072, 2.629628906250007, 2.7296289062500065, 2.829628906250006, 2.929628906250006, 3.0296289062500055, 3.129628906250005, 3.2296289062500048, 3.3296289062500044, 3.429628906250004, 3.5296289062500037, 3.6296289062500033, 3.729628906250003, 3.8296289062500026, 3.9296289062500023, 4.029628906250002, 4.129628906250002, 4.229628906250001, 4.329628906250001, 4.4296289062500005, 5.72],
NON_UNIFRM = False, QUANT = False, ADJ_BIN = False,
LABELS = True, EXTEND = True,
title=r'$m_c=%s~b=%s$'%(mc,bval), #'Yuha Desert',
save = '%s/YuhaDeserts_prod'%path_save, dpi = 150) #67)
ax.set_title('mc=%s,~bval=%s'%(mc,bval),fontsize='xx-small')
#print mlist
#np.save('%s/alpha'%DIR_OUTPT_figs,alpha_p)
In [38]:
def GetMlist( nij_trig, swarm_lohi ):
#--- remove pairs with r > cs.t
df_trig = nij_trig.copy()
#--- parents' magnitudes
m_list = GetParentsMagnitudes( swarm_lohi, df_trig )
#--- number of daughters for each mother
df_nn = df_trig.groupby(by='Parent_id').count()['Event_id']
#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
# print swarm_lohi['magnitude'].loc[df_zeros.index]
# print df_zeros
#--- parents' magnitudes including those with zero child
mm = swarm_lohi['magnitude'].loc[df_zeros.index].tolist() + \
swarm_lohi['magnitude'].loc[df_nn.index].tolist() #swarm_lohi[swarm_lohi['magnitude']>=mc]['magnitude']
mm = np.array(mm)
# mm=list(set(mm))
mm.sort()
mlist_with_zero=np.array(mm)
return mlist_with_zero
def GetMhi10( mlist ):
nmax = 10
N=len(mlist)
hist, edges = np.histogram(mlist,bins=1024*10,density=None)
assert np.sum(hist) == N
# plt.plot(edges,N-np.append(np.cumsum(hist),[np.cumsum(hist)[-1]]),drawstyle='steps-post')
# plt.plot([2.939517578125,2.939517578125],[0,10],'-')
# plt.plot([2.93984375,2.93984375],[0,10],'-')
# plt.ylim(0,11)
# plt.xlim(2.75,4)
# plt.show()
# Quantiles = np.concatenate([(1.0*np.arange(nn)/nn)[1:],np.array([1])])
q = 1.0 - 1.0 * nmax / N
# print Quantiles
# m_q=[np.min(mlist)-1e-6]
# for q in Quantiles:
# print N - np.cumsum(hist)
return edges[:-1][np.cumsum(hist)<=q*N][-1], edges[1:][np.cumsum(hist)<=q*N][-1]
mlist_with_zero = GetMlist( nij_trig, swarm_lohi )
mhi = np.mean( GetMhi10( mlist_with_zero ) )
#print len(mlist_with_zero[mlist_with_zero>mhi])
m_max = mbins[-1]
a = np.append(np.array([m_max]),np.arange(mhi,mc-0.1,-0.1))
a.sort()
print list(a)
In [39]:
# n=1024
# a = np.random.normal(size=n) #np.random.random(n)
# bin_edges = stats.bayesian_blocks(a) #,fitness='measures') #,dt=0.1)
# # bin_edges=np.histogram_bin_edges(a, bins='auto')
# nbin=len(bin_edges)-1
In [40]:
# plt.scatter(a,np.zeros(n),marker='x',alpha=0.1)
# plt.plot(bin_edges,np.concatenate([np.arange(nbin),[nbin-1]],axis=0),drawstyle='steps-post')
# plt.show()
In [41]:
n = len(mbins[mbins>3.0]) #16
MHI = mbins[mbins>3.0] #np.linspace(3,6,n)
pval_lin = np.zeros(n)
pval_nonlin = np.zeros(n)
for mhi, index in zip(MHI,xrange(sys.maxint)):
pval_lin[index], pval_nonlin[index],junk = ProdRelation( nij_trig, swarm_lohi,
NON_UNIFRM = True,
mlo=mc-0.05, mhi=mhi+0.05 )
In [42]:
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.set_yscale('log')
ax.plot(MHI,pval_lin,'-o',label='linear')
ax.plot(MHI,pval_nonlin,'-s',label='nonlinear')
ax.legend(fontsize='xx-small',frameon=False)
ax.set_xlabel(r'$m_{max}$')
ax.set_ylabel(r'$p$')
plt.savefig('%s/p.png'%DIR_OUTPT_figs,dpi=100,bbox_inches='tight')
plt.show()
In [43]:
pval_nonlin
In [44]:
MHI
multiple realizations¶
In [45]:
from IPython.display import Image
alpha_arr=[]
for irun in xrange(16):
path = 'runYuha/Run%s'%irun
a=np.load('%s/alpha.npy'%path)
alpha_arr.append(a.item())
#--- load image
display(Image(filename='%s/prod.png'%path))
df=pd.DataFrame({'mean':[np.mean(alpha_arr)],'std_dev':np.std(alpha_arr),
'std_err':[np.std(alpha_arr)/len(alpha_arr)**.5]},
index=[r'$\alpha$'])
display(df)
df['mean']+2*df['std_dev'],df['mean']-2*df['std_dev']
In [46]:
N_as
Out[46]:
In [47]:
#--- choose bins for rho and lambda
#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi.reset_index(inplace=True,drop=True)
#--- remove pairs with r > cs.t
#cs = 3 * (24*3600) #--- km per day
df_trig = nij_trig.copy()
#--- parents' magnitudes
m_list = GetParentsMagnitudes( swarm_lohi, df_trig )
#--- number of daughters for each mother
df_nn = df_trig.groupby(by='Parent_id').count()['Event_id']
#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
#--- scattered plot
fig = plt.figure(figsize=(4,4),frameon=True)
ax = fig.add_subplot(111)
#ax.set_xlim(mc,swarm_lohi['magnitude'].max()+0.5)
m_min = min(m_list)-1.0e-6
m_max = max(m_list)+1.0e-6 #--- exclude the largest event
ax.set_xlim(m_min,np.ceil(m_max))
#ax.set_ylim(0.5,1000)
ax.plot( swarm_lohi['magnitude'].loc[df_nn.index],
df_nn,
'o',color='black',alpha=.1)
ax.plot( swarm_lohi['magnitude'].loc[df_zeros.index],
df_zeros,
'x',color='gray')
plt.yscale('log')
#plt.xlabel('$m$')
#plt.ylabel(r'$N_{as}$')
#--- choose m0<m<m1
dm = 0.7
m_min = mc-1.0e-6
m_max = 5.0 #swarm_lohi['magnitude'].max()+1.0e-6 #--- exclud the largest one
n = int( ( m_max - m_min ) / dm )
m = list(np.linspace(m_min,m_max,n))
#df_parents = df_trig.groupby(by='Parent_id').groups #--- parent events
#ds = swarm_lohi['magnitude'].loc[df_parents.keys()] #--- mag. of parent events
#m=AdjustBins(ds,m)
#m.pop(-2)
m+=[6.0]
Mlist=m
m_mean=np.array(np.zeros(len(m) - 1))
m_err=np.array(np.zeros(len(m) - 1))
N_as = np.array(np.zeros(len(m) - 1))
err_N = np.array(np.zeros(len(m) - 1))
for i in xrange( len(m) - 1 ):
try:
m0 = m[ i ]
m1 = m[ i + 1 ]
df, m_mean[ i ], m_err[i] = GetPairsWithSpecifiedParentMag((m0,m1),
swarm_lohi,df_trig) #--- get parent with m0<m<m1
df.reset_index( inplace = True, drop = True )
# pdb.set_trace()
#--- number of daughters for each mother
df_n = df.groupby(by='Parent_id').count()['Event_id']
#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
subset = swarm_lohi.loc[list_of_zero_afs]
list_of_zero_afs = subset[(subset['magnitude'] >= m0) & (subset['magnitude'] <= m1)].index
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
#--- mean
# print [m0,m1],list_of_zero_afs
# pdb.set_trace()
N_as[ i ] = 1.0*(df_n.sum()+ df_zeros.sum())/(len(df_n)+len(df_zeros))
err_N[ i ] = (N_as[ i ]/(len(df_n)+len(df_zeros)))**0.5 #--- poisson
# err_N[ i ]=np.std(df_n.tolist()+df_zeros.tolist())/(len(df_n)+len(df_zeros))**0.5
except:
# traceback.print_exc()
continue
#ax.plot(m, np.append(N_as,N_as[-1]), linestyle='-', drawstyle='steps-post',color='red')
ax.errorbar(m_mean,N_as,
yerr=err_N,xerr=m_err,
markerfacecolor='white',
fmt='o', markersize=8,color='red',markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1)
#--- add major xticks
xmin=np.ceil(ax.axis()[0])
xmax=np.floor(ax.axis()[1])
nbin = xmax - xmin
ax.set_xticks(np.linspace(xmin,xmax,int(nbin)+1))
DrawFrame(ax, (0.15,0.08),(0.14,0.06),0.1,LOG_Y=True)
try:
from scipy.optimize import curve_fit
def func(x,c,alpha):
return c*10**(alpha*x)
popt, pconv = curve_fit(func,m_mean[m_mean>=mc],N_as[m_mean>=mc])
mm = swarm_lohi['magnitude'].loc[df_zeros.index].tolist() + \
swarm_lohi['magnitude'].loc[df_nn.index].tolist() #swarm_lohi[swarm_lohi['magnitude']>=mc]['magnitude']
mm = np.array(mm)
mm=list(set(mm))
mm.sort()
mm=np.array(mm)
ax.plot(mm,func(mm,*popt),'-.g')
print popt
print pconv**0.5
except:
pass
#plt.savefig('%s/prod.png'%DIR_OUTPT_figs,dpi=150,bbox_inches='tight')
fig = plt.figure(figsize=(4,4),frameon=True)
ax = fig.add_subplot(111)
ax.plot(m, np.append(N_as,N_as[-1]), linestyle='-', drawstyle='steps-post',color='red')
ax.errorbar(m_mean,N_as,
yerr=err_N,xerr=m_err,
markerfacecolor='white',
fmt='o', markersize=8,color='red',markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1)
ax.set_xlim(m_min,np.ceil(m_max))
#ax.set_ylim(0.5,1000)
ax.set_yscale('log')
ax.set_xticks(np.linspace(xmin,xmax,int(nbin)+1))
plt.show()
In [48]:
#--- density plots: rho(r)
#--- split based on the mainshock's magnitude
import numpy as np
import math
def AdjustBins(data,bins):
hist, xedges=np.histogram(data, #--- histogram
bins=bins)
tmp = [xedges[0]]
for i,j in zip(xedges[1:],hist): #--- expand bins
if j < 2:
continue
tmp.append(i)
# tmp.append(xedges[-1])
if hist[-1] < 2:
tmp.pop(-2)
return tmp
# def GetPairsWithSpecifiedParentMag((m0,m1),catalog,df):
# df_parents = df.groupby(by='Parent_id').groups #--- parent events
# ds = catalog['magnitude'].loc[df_parents.keys()] #--- mag. of parent events
# ds_m=ds[(m0 <= ds) & (ds < m1)] #--- parent events with m0<m<m1
# df_m=df[pd.DataFrame([df['Parent_id'] ==k for k in ds_m.index]).any()] #--- data frame
# return df_m, np.mean( ds_m )
# def AddDist( df_trig, df_complete ):
# x=df_complete.loc[ df_trig['Event_id'] ]['r(m)']
# y=df_complete.loc[ df_trig['Parent_id'] ]['r(m)']
# df_trig['R'] = pd.Series(numpy.linalg.norm(x-y)))
# assert len ( df_trig[ df_trig['R'] == 0.0 ] ) == 0, '%s'%display( df_trig[ df_trig['R'] == 0.0 ] )
# def AddTime( df_trig, df_complete ):
# prefact = 1.0 / ( 24.0 * 60 * 60 ) #--- daily
# x=df_complete.loc[ df_trig['Event_id'] ]
# y=df_complete.loc[ df_trig['Parent_id'] ]
# x.reset_index(inplace=True)
# y.reset_index(inplace=True)
# df_trig['T'] = x['date']-y['date']
# df_trig['T'] = df_trig['T'].apply(lambda x: x.total_seconds() * prefact)
# def RemovePair( df, cs ):
# return df [ df[ 'R' ] <= df[ 'T' ] * cs ]
#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi.reset_index(inplace=True,drop=True)
#--- choose m0<m<m1
X={}
rho_dict={}
m_range={}
Err={}
m_mean={}
for i in xrange( len(m) - 1 ):
try:
m0 = m[ i ]
m1 = m[ i + 1 ]
df_trig, m_mean[ i ],junk = GetPairsWithSpecifiedParentMag((m0,m1),swarm_lohi,nij_trig) #--- get parent with m0<m<m1
df_trig.reset_index( inplace = True, drop = True )
#--- rho plots
n_decades = int( math.ceil( log(df_trig['R'].max()/df_trig['R'].min(),10) ) )
nbin_per_decade = 4
nbins = nbin_per_decade * n_decades
# if i == 3: #--- remove
df_trig=df_trig[df_trig['R'] <= 7.0] #--- r < 10
bins = np.logspace(log(df_trig['R'].min(),10)-1.e-10,
log(df_trig['R'].max(),10)+1.e-10,nbins)
bins = AdjustBins(df_trig['R'],bins)
rho, xedges=np.histogram(df_trig['R'],
density=True,
bins=bins)
hist, xedges=np.histogram(df_trig['R'],
bins=bins)
x_bins, junk=np.histogram(df_trig['R'], #--- histogram
bins=bins,
weights = df_trig['R'])
#--- rm empty bins
x_bins = x_bins[ hist != 0 ]
rho = rho[ hist != 0 ]
hist = hist[ hist != 0 ]
#--- poisson err
err = rho/np.sqrt(hist)
x_bins /= hist
#--- scattered plot
plt.figure(figsize=(8,4))
plt.subplot(1,2,1)
plt.title('$%2.1f<m<%2.1f$'%(m0,m1))
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$t$')
plt.ylabel('$r$')
dt_min=nij_trig['T'].min()
dt_max=nij_trig['T'].max()
dr_min=nij_trig['R'].min()
dr_max=nij_trig['R'].max()
plt.ylim(dr_min/2,dr_max*2)
plt.xlim(dt_min/2,dt_max*2)
plt.scatter(df_trig['T'],df_trig['R'],alpha=0.1)
plt.plot([dt_min,dt_min],[dr_min,dr_max],
'-',color='blue')
plt.plot([dt_max,dt_max],[dr_min,dr_max],
'-',color='blue')
plt.plot([dt_min,dt_max],[dr_min,dr_min],
'-',color='blue')
plt.plot([dt_min,dt_max],[dr_max,dr_max],
'-',color='blue')
plt.plot([dt_min,dt_max],[dt_min*cs,dt_max*cs],
'-',color='red')
#--- rho
plt.subplot(1,2,2)
# plt.xlim(1e-3,1e2)
# plt.ylim(1e-3,1e2)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$r$')
plt.ylabel('$P(r)$')
plt.errorbar(x_bins,rho, yerr=err)
X[i] = x_bins[:]
rho_dict[i]=rho[:]
m_range[i]=np.array([m0,m1])
Err[i]=err[:]
except:
traceback.print_exc()
continue
plt.show()
$r_{rup}$¶
In [49]:
#--- density plots: rho(r)
#--- split based on the mainshock's magnitude
import numpy as np
import math
#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi.reset_index(inplace=True,drop=True)
XX={}
rho_dictt={}
Errr={}
m_meann={}
r_rupp = {}
groups_dic=nij_trig.groupby(by='Parent_id').groups
for parent_id in groups_dic:
df_trig = nij_trig.loc[groups_dic[parent_id]]
df_trig.reset_index( inplace = True, drop = True )
if len(df_trig) < 2:
continue
i = parent_id
if np.std(df_trig['R']) == 0.0:
continue
m_meann[i]=swarm_lohi['magnitude'].loc[i]
r_rupp[i]=np.std(df_trig['R'])
m_meann=pd.Series(m_meann)
r_rupp=pd.Series(r_rupp)
#--- fit
def func(x,alpha,c):
return c*10**(alpha*x)
popt = np.polyfit(m_meann,np.log10(r_rupp), 1)
popt[1] = 10**(popt[1])
#--- plot
fig = plt.figure(figsize=(4,4))
ax=fig.add_subplot(111)
#ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('$m$')
ax.set_ylabel('$r_{rup}$[km]')
#ax.set_xlim(2,6)
#ax.set_ylim(1e-4,1e2)
ax.scatter(m_meann,r_rupp,alpha=0.2)
ax.plot(m_meann,func(m_meann,*popt),'-.r',zorder=3,
label='$r_{rup}\propto 10^{%3.2f m},r_c=%2.1e$'%(popt[0],popt[1]))
ax.legend(fontsize='xx-small',frameon=False)
#ax.set_title('q=%d%%,df=%s'%(quantile*100,Df),fontsize='xx-small')
sigma=popt[0]
rc=popt[1] #10**np.ceil(np.log10(popt[1]))
##plt.savefig('%s/rupture_mag.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight',transparent=True)
plt.show()
$r_{rup}$ by fitting: run twice¶
In [50]:
from scipy.optimize import curve_fit
#def func(x,pref, xmean,sigma):
# return pref * np.exp(-((x-xmean)/sigma)**2)
def func(x,a,b1,b2,rc):
return a*x**b1/(1.0+(x/rc)**b2)
#def func2(x,c,alpha):
# return c*10**(alpha*x)
#--- fit
def func2(x,alpha,c):
return c*10**(alpha*x)
colors = ['black','red','green','blue','cyan','brown','violet','magenta','orange','yellow']
fig=plt.figure(figsize=(8,4))
ax=fig.add_subplot(121)
#ax.axis([2e-3,2e-3*1e1/1e-3,1e-3,1e1])
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('$r$')
ax.set_ylabel('$P(r)$')
r_peak=[]
mm_mean=[]
p0 = (1.0,Df-1,2.0,rc*10**(sigma*m_mean[0])) #--- initialize
for key, col in zip(X,colors):
try:
# if key > 6:
# continue
# print np.mean(X[key])
ax.errorbar(X[key],rho_dict[key],
yerr=Err[key],
linestyle='None', color=col,
fmt='-o', markersize=8
)
popt, pcov = curve_fit(func, X[key], rho_dict[key],
sigma=Err[key],
# bounds=((0.0,0.0,0.0), (np.inf,np.inf,np.inf)),
p0=p0
)
mm_mean+=[m_mean[key]]
p0 = popt[:]
(a,b1,b2,rc)=popt
r_peak += [rc]
ax.plot(X[key], func(X[key], *popt),'-',color=col)#,label='fit: r=%3.2f' % (popt[1]))
ax.legend()
except:
traceback.print_exc()
# print key
continue
try:
ax=fig.add_subplot(122)
ax.set_yscale('log')
ax.set_xlabel('$m$')
ax.set_ylabel('$r_{rup}$')
ax.plot(mm_mean,r_peak,'o')
ylo=10**np.floor(np.log10(np.min(r_peak)))
yhi=10**np.ceil(np.log10(np.max(r_peak)))
ax.set_ylim(ylo,yhi)
#
popt = np.polyfit( mm_mean,np.log10(r_peak),1 )
sigma=popt[0]
rc=10**popt[1] #10**np.ceil(np.log10(popt[1]))
rc0 = rc
#
ax.plot(mm_mean,func2( np.array(mm_mean),sigma,rc),'--r',
label='$r_{rup}\propto 10^{%3.2f m},r_c=%2.1e$'%(sigma,rc))
ax.legend(fontsize='xx-small',frameon=False)
fig.tight_layout()
except:
traceback.print_exc()
pass
plt.show()
fig.savefig('%s/r_rupt.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight')
np.save('%s/sigma'%DIR_OUTPT_figs,sigma)
plt.show()
different realizations¶
In [51]:
alpha_arr=[]
for irun in xrange(16):
path = 'runYuha/Run%s'%irun
a=np.load('%s/sigma.npy'%path)
alpha_arr.append(a.item())
#--- load image
display(Image(filename='%s/r_rupt.png'%path))
df=pd.DataFrame({'mean':[np.mean(alpha_arr)], 'std_dev':np.std(alpha_arr),
'std_err':[np.std(alpha_arr)/len(alpha_arr)**.5]},
index=[r'$\sigma$'])
display(df)
df['mean']+2*df['std_dev'],df['mean']-2*df['std_dev']
In [52]:
#--- density plots: rho(r) 2nd
plt.figure(figsize=(15,5))
#--- plot rho
plt.subplot(1,3,1)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'$r$')
plt.ylabel(r'$P(r)$')
plt.xlim(1e-3,1e2)
plt.ylim(1e-4,1e1)
for key in X:
# if key < 3:
# continue
plt.errorbar(X[key],rho_dict[key],
yerr=Err[key],
fmt='-o', markersize=8
)
#--- rescaled in x and y
plt.subplot(1,3,2)
alpha=sigma
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'$R=r/10^{%3.2f m}$'%sigma)
plt.ylabel(r'$P(R)=P(r)\times 10^{%3.2f m}$'%sigma)
plt.xlim(1e-5,1e0)
plt.ylim(1e-2,1e3)
xdata = []
ydata = []
y_err=[]
for key in X:
if key > 3:
continue
xdata += list(X[key]/10**(sigma*np.mean(m_range[key])))
ydata += list(rho_dict[key]*10**(sigma*np.mean(m_range[key])))
y_err += list(Err[key]*10**(alpha*m_mean[key]))
plt.errorbar(X[key]/10**(sigma*m_mean[key]),
rho_dict[key]*10**(sigma*m_mean[key]),
yerr=Err[key]*10**(sigma*m_mean[key]),
label='%2.1f-%2.1f'%(m_range[key][0],m_range[key][1]),
fmt='-o', markersize=8
)
plt.legend(title='m',fontsize='small',loc='upper right',frameon=False,bbox_to_anchor=(0.44, .64))
#--- rescaled by r^-\nu
nu = 1.7
plt.subplot(133)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'$R$')
plt.ylabel(r'$P(R)\times R^{%s}$'%nu)
#plt.ylabel(r'$\lambda(t) \times t^p$')
plt.xlim(1e-4,1e0)
#plt.ylim(1e-3,1e5)
KEYS=X.keys()
KEYS.sort()
for key, index in zip(KEYS,xrange(100)):
# if key < 3:
# continue
x_new = X[key]/10**(sigma*m_mean[key])
y_new = rho_dict[key]*10**(sigma*m_mean[key])
y2_new = Err[key]*10**(sigma*m_mean[key])
plt.errorbar(x_new,
y_new*x_new**nu,
yerr=y2_new*x_new**nu,
fmt='-o', markersize=8
)
plt.tight_layout()
#plt.savefig('%s/rho_r.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight')
plt.show()
power-law fit¶
$\rho(r)\propto \frac{r^{\nu_1}}{1+(r/r_c)^{\nu_2}}$¶
In [53]:
#--- fitting
from scipy.optimize import curve_fit
def func(x,a,b1,b2,rc):
return a*x**b1/(1.0+(x/rc)**b2)
xdata = np.array(xdata)
xxdata = xdata[:] #[xdata>rc]
#
ydata = np.array(ydata)
yydata = ydata[:] #[xdata>rc]
#
Y_err=np.array(y_err)
#
ss=zip(xxdata,yydata,Y_err)
ss.sort()
xxdata=[i[0]for i in ss]
yydata=[i[1]for i in ss]
Y_err = [i[2] for i in ss]
popt, pcov = curve_fit(func, xxdata, yydata,
p0=(1.0,Df-1,2.0,rc),sigma=Y_err) #, bounds=((0.0,0.0), (np.inf,10.0)))
(a,b1,b2,rc)=popt
nu = -b2
rho_pref = func(rc,a,b1,b2,rc)
#--- plot
fig = plt.figure(figsize=(4,4))
ax=fig.add_subplot(111)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('$r[km]$')
ax.set_ylabel(r'$\rho(r)$')
#ax.set_xlim(2,6)
#ax.set_ylim(1e-4,1e2)
ax.scatter(xxdata,yydata)
# #--- fit
ax.plot(xxdata,func(xxdata,*popt),
'-.r',zorder=3,
label=r'$\rho(r)\propto \frac{r^{%3.2f}}{1+(r/r_c)^{%3.2f}},r_c={%2.1e}$'%(b1,b2,rc))
ax.legend(fontsize='xx-small',frameon=False)
#ax.set_title('q=%d%%,df=%s'%(quantile*100,Df),fontsize='xx-small')
###plt.savefig('%s/rupture_mag.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight',transparent=True)
#np.save('%s/nu'%DIR_OUTPT_figs,nu+b1)
plt.show()
In [54]:
from scipy.stats import chi2
#--- fitting
from scipy.optimize import curve_fit
def func(x,b,a):
return a*x**b
xdata = np.array(xdata)
xxdata = xdata[xdata>rc]
ydata = np.array(ydata)
yydata = ydata[xdata>rc]
Y_err=np.array(y_err)
Y_err = Y_err[xdata>rc]
#popt, pcov = curve_fit(func, xxdata, yydata, p0=(1.0,1.0), bounds=((0.0,0.0), (np.inf,10.0)))
popt = np.polyfit(np.log10(xxdata),np.log10(yydata), 1)
popt[1] = 10**(popt[1])
nu = popt[0]
rho_pref = popt[1]
#--- goodness of fit
n=len(xxdata)
m=2 #--- number of coefficients
xsq = np.sum(((yydata-func(xxdata,*popt))/Y_err)**2)
pval=1-chi2.cdf(xsq,n-m)
#print pval
#--- plot
fig = plt.figure(figsize=(4,4))
ax=fig.add_subplot(111)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('$r[km]$')
ax.set_ylabel(r'$\rho(r)$')
#ax.set_xlim(2,6)
#ax.set_ylim(1e-4,1e2)
ax.scatter(xdata,ydata)
#--- fit
ax.plot(xxdata,func(xxdata,*popt),
'-.r',zorder=3,
label=r'$\rho(r)\propto r^{%3.2f}$'%popt[0])
ax.legend(fontsize='xx-small',frameon=False)
#ax.set_title('q=%d%%,df=%s'%(quantile*100,Df),fontsize='xx-small')
##plt.savefig('%s/rupture_mag.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight',transparent=True)
plt.show()
new goodness of fit for yuha¶
In [55]:
from scipy.stats import chi2
#--- fitting
from scipy.optimize import curve_fit
def func(x,a,b1,b2,rc):
return a*x**b1/(1.0+(x/rc)**b2)
def func2(x,b,a):
return a*x**b
pval_max = 0.0
for sigmaa in np.linspace(0.2,0.5,32):
#--- rescale
x = []
y = []
ye = []
for key in X:
if key > 3:
continue
x += list(X[key]/10**(sigmaa*np.mean(m_range[key])))
y += list(rho_dict[key]*10**(sigmaa*np.mean(m_range[key])))
ye += list(Err[key]*10**(sigmaa*m_mean[key]))
#--- initialize
x = np.array(x)
y = np.array(y)
ye=np.array(ye)
#
ss=zip(x,y,ye)
ss.sort()
x=np.array([i[0]for i in ss])
y=np.array([i[1]for i in ss])
ye = np.array([i[2] for i in ss])
#--- fit
popt0, pcov = curve_fit(func, x, y,
p0=(1.0,Df-1,2.0,rc0),sigma=ye) #, bounds=((0.0,0.0), (np.inf,10.0)))
(a,b1,b2,rc)=popt0
# print rc
#--- goodness of fit
# print rc,x
xxdata = x[x>rc]
yydata = y[x>rc]
Y_err = ye[x>rc]
popt1 = np.polyfit(np.log10(xxdata),np.log10(yydata), 1)
popt1[1] = 10**(popt1[1])
nu = popt1[0]
rho_pref = popt1[1]
n=len(xxdata)
m=2 #--- number of coefficients
xsq = np.sum(((yydata-func2(xxdata,*popt1))/Y_err)**2)
pval=1-chi2.cdf(xsq,n-m)
# print pval
if pval > pval_max:
pval_max = pval
Xdata = x[:]
Ydata = y[:]
popt = popt1[:]
sigma = sigmaa
nu = popt1[0]
rho_pref = popt1[1]
RC = rc
rc = RC
#--- plot
fig = plt.figure(figsize=(4,4))
ax=fig.add_subplot(111)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('$r[km]$')
ax.set_ylabel(r'$\rho(r)$')
#ax.set_xlim(2,6)
#ax.set_ylim(1e-4,1e2)
ax.scatter(Xdata,Ydata)
# #--- fit
ax.plot(Xdata[Xdata>RC],func2(Xdata[Xdata>RC],*popt),
'-.r',zorder=3,
label=r'$\rho(r)\propto r^{%3.2f}$'%nu)
ax.legend(fontsize='xx-small',frameon=False)
ax.set_title(r'$q=%d,df=%s,\sigma=%3.2f,p=%3.2e$'%(quantile*100,Df,sigma,pval_max),fontsize='xx-small')
#plt.savefig('%s/fit.png'%DIR_OUTPT_figs,dpi=150,bbox_inches='tight',transparent=True)
plt.show()
np.save('%s/nu'%DIR_OUTPT_figs,nu)
np.save('%s/sigma_max'%DIR_OUTPT_figs,sigma)
np.save('%s/pval'%DIR_OUTPT_figs,pval_max)
different realizations¶
In [56]:
import numpy as np
import pandas as pd
alpha_arr=[]
PVAL = []
SIGMA = []
for irun in xrange(16):
a=np.load('runYuha/Run%s/nu.npy'%(irun))
b=np.load('runYuha/Run%s/pval.npy'%(irun))
c=np.load('runYuha/Run%s/sigma_max.npy'%(irun))
# print irun,a.item()
# if b.item() < 0.001: #--- remove p < 0.001
# continue
alpha_arr.append(a.item())
PVAL.append([b.item(),irun])
SIGMA.append(c.item())
df=pd.DataFrame({'mean':[np.mean(alpha_arr),np.mean(SIGMA)],
'std_dev':[np.std(alpha_arr),np.std(SIGMA)],
'std_err':[np.std(alpha_arr)/len(alpha_arr)**.5,np.std(SIGMA)/len(SIGMA)**.5]},
index=[r'$\nu$',r'$\sigma$'])
display(df)
print r'%2.1f<\nu<%2.1f'%((df['mean']+2*df['std_dev']).iloc[0],(df['mean']-2*df['std_dev']).iloc[0])
print r'%2.1f<\sigma<%2.1f'%((df['mean']+2*df['std_dev']).iloc[1],(df['mean']-2*df['std_dev']).iloc[1])
In [57]:
PVAL.sort(reverse=True)
PVAL
Out[57]:
In [58]:
#--- density plots: rho(r) (paper version)
from matplotlib import rc as rcc
import matplotlib as mpl
mpl.rcParams.update(mpl.rcParamsDefault)
import matplotlib.ticker
rcc('text', usetex=True)
font = {'size' : 20}
matplotlib.rc('font', **font)
def returnCMtoInch(x,y):
return (x/2.54,y/2.54)
fig, ax = plt.subplots(figsize=(4,4))
#sigma=alpha
xdata = []
ydata = []
colors = ['black','red','green','blue','cyan','brown','violet','magenta','orange','yellow']
fillstyles=['white',None,'white',None,'white',None,'white',None,'white',None,'white',None,'white',None,'white',None]
markers=['o','s','D','^','<','>','v']
assert len(colors) >= len(X)
for key,col,mk,fl in zip(X,colors,markers,fillstyles):
xdata += list(X[key]/10**(sigma*np.mean(m_range[key])))
ydata += list(rho_dict[key]*10**(sigma*np.mean(m_range[key])))
ax.errorbar(X[key]/10**(sigma*m_mean[key]),
rho_dict[key]*10**(sigma*m_mean[key]),
yerr=Err[key]*10**(sigma*m_mean[key]),
label='%2.1f-%2.1f'%(m_range[key][0],m_range[key][1]),
fmt='-o', markersize=8,color=col,marker=mk,markerfacecolor=fl,markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1
)
(ylo,yhi)=(10**np.floor(np.log10(np.min(ydata))),\
10**np.ceil(np.log10(np.max(ydata))))
xhi=10**np.ceil(np.log10(np.max(xdata)))
ax.axis([xhi/(yhi/ylo),xhi,ylo,yhi])
ax.loglog()
#--- add major xticks
xmin=np.ceil(np.log10(ax.axis()[0]))
xmax=np.floor(np.log10(ax.axis()[1]))
nbin = xmax - xmin
ax.set_xticks(np.logspace(xmin,xmax,int(nbin)+1))
#--- put minor bins
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
#plt.xlabel(r'$R=r/10^{%s m}$'%alpha)
#plt.ylabel(r'$P(R)=P(r)\times 10^{%s m}$'%alpha)
#fig.legend(title='m',prop=fontP,loc='upper right',frameon=False)#, bbox_to_anchor=(0.4, .6))
ax.tick_params(axis='y',left=True, right=True,which='both')
ax.tick_params(axis='x',bottom=True, top=True,which='both')
#--- plot power-law
(xlo,xhi)=(rc,100*rc)
yhi=rho_pref*rc**nu/10**.5
ax.plot([xlo,xhi],[yhi, yhi* (xhi/xlo)**(nu)],'r-.',
label=r'$\rho(r)\propto r^{%3.2f},r_{rup}\propto 10^{%3.2f m}$'%(nu,sigma),
linewidth=1)
#--- rescaled
#ax.scatter(xdata,1e2*np.array(ydata)/np.array(xdata)**(nu))
ax.legend(fontsize='x-small',frameon=False)
ax.set_xlabel(r'$r/r_{rup}$')
#ax.set_ylabel(r'$\rho(r)\times r_{rup}$')
ax.set_title('mc=%s,~b=%s'%(mc,bval),fontsize='xx-small')
#DrawFrame(ax, (0.22,0.06),(0.14,0.06),0.1,LOG_X=True,LOG_Y=True)
#plt.tight_layout()
plt.savefig('%s/rho_rescaled.png'%DIR_OUTPT_figs,dpi=75*2,bbox_inches='tight')
plt.show()
different realizations¶
In [59]:
from IPython.display import Image
for irun in xrange(16):
path = 'runYuha/Run%s'%irun
#--- load image
display(Image(filename='%s/rho_rescaled.png'%path))
In [60]:
#--- density plots: rho(r) (paper version)
from matplotlib import rc as rcc
import matplotlib as mpl
mpl.rcParams.update(mpl.rcParamsDefault)
import matplotlib.ticker
rcc('text', usetex=True)
font = {'size' : 20}
matplotlib.rc('font', **font)
#plt.style.use('seaborn')
#plt.rcParams.update({'lines.markeredgewidth': 1})
def returnCMtoInch(x,y):
return (x/2.54,y/2.54)
fig, ax = plt.subplots(figsize=(4,4))
#sigma=alpha
xdata = []
ydata = []
y_err = []
colors = ['black','red','green','blue','cyan','brown','grey','magenta','orange','yellow']
fillstyles=['white',None,'white',None,'white',None,'white',None,'white',None,'white',None,'white',None,'white',None]
markers=['o','s','D','^','<','>','v']
assert len(colors) >= len(X)
for key,col,mk,fl in zip(X,colors,markers,fillstyles):
xdata += list(X[key]/10**(sigma*np.mean(m_range[key])))
ydata += list(rho_dict[key]*10**(sigma*np.mean(m_range[key])))
y_err += list(Err[key]*10**(sigma*m_mean[key]))
ax.errorbar(X[key]/10**(sigma*m_mean[key]),
rho_dict[key]*10**(sigma*m_mean[key]),
yerr=Err[key]*10**(sigma*m_mean[key]),
# label='%2.1f-%2.1f'%(m_range[key][0],m_range[key][1]),
fmt='-o', markersize=8,color=col,marker=mk,markerfacecolor=fl,markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1
)
xdata=np.array(xdata)
ydata=np.array(ydata)
y_err=np.array(y_err)
(ylo,yhi)=(10**np.floor(np.log10(np.min(ydata-y_err))),\
# 10**np.ceil(np.log10(np.max(ydata))))
np.max(ydata))
#(xlo,xhi)=(10**np.floor(np.log10(np.min(xdata))),\
# 10**np.ceil(np.log10(np.max(xdata))))
xc=(np.min(xdata)*np.max(xdata))**.5
xlo = xc/(yhi/ylo)**.5
xhi = xc*(yhi/ylo)**.5
ax.axis([xlo,xhi,ylo,yhi])
ax.loglog()
#--- add major xticks
xmin=np.ceil(np.log10(ax.axis()[0]))
xmax=np.floor(np.log10(ax.axis()[1])) #-1
nbin = xmax - xmin
ax.set_xticks(np.logspace(xmin,xmax,int(nbin)+1))
#--- put minor bins
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
#plt.xlabel(r'$R=r/10^{%s m}$'%alpha)
#plt.ylabel(r'$P(R)=P(r)\times 10^{%s m}$'%alpha)
#fig.legend(title='m',prop=fontP,loc='upper right',frameon=False)#, bbox_to_anchor=(0.4, .6))
ax.tick_params(axis='y',left=True, right=True,which='both')
ax.tick_params(axis='x',bottom=True, top=True,which='both')
#--- plot power-law
#(xlo,xhi)=(rc,100*rc)
#yhi=rho_pref/10**.5
#ax.plot([xlo,xhi],[yhi, yhi* (xhi/xlo)**(nu)],'r-.',
# label=r'$\rho(r)\propto r^{%3.2f},r_{rup}\propto 10^{%3.2f m}$'%(nu+b1,sigma),
# linewidth=1)
#--- plot power-law
(xlo,xhi)=(rc,100*rc)
yhi=rho_pref*rc**nu/10**.5
ax.plot([xlo,xhi],[yhi, yhi* (xhi/xlo)**(nu)],'r-.',
# label=r'$\rho(r)\propto r^{%3.2f},r_{rup}\propto 10^{%3.2f m}$'%(nu,sigma),
linewidth=1)
#ax.legend(fontsize='x-small',frameon=False)
#ax.set_xlabel(r'$r/r_{rup}$')
#ax.set_ylabel(r'$\rho(r)\times r_{rup}$')
#ax.set_title('q=%d%%,df=%s'%(quantile*100,Df),fontsize='xx-small')
DrawFrame(ax, (0.22,0.06),(0.14,0.06),0.1,LOG_X=True,LOG_Y=True)
#plt.tight_layout()
plt.savefig('%s/rho_r_2nd.png'%DIR_OUTPT_figs,dpi=300,bbox_inches='tight')
plt.show()
In [61]:
#--- density plots: rho(r) (paper version 2nd)
#plt.style.use('seaborn')
#plt.rcParams.update({'lines.markeredgewidth': 1})
mpl.rcParams.update(mpl.rcParamsDefault)
rcc('text', usetex=True)
font = {'size' : 20}
matplotlib.rc('font', **font)
def returnCMtoInch(x,y):
return (x/2.54,y/2.54)
fig, ax = plt.subplots(figsize=(4,4))
#---
ax.spines['right'].set_visible(False) #--- half open
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
#---
#ax.axis([1e-3,1e2,1e-4,1e1])
ax.loglog()
#ax.set_xticks(np.logspace(-1,3,4))
#plt.xlabel(r'$R=r/10^{%s m}$'%alpha)
#plt.ylabel(r'$P(R)=P(r)\times 10^{%s m}$'%alpha)
#plt.xlim(1e-5,1e0)
#plt.ylim(1e-5,1e0)
#colors = ['black','red','green','blue','cyan','brown']
#fillstyles=['white',None,'white',None,'white',None,'white']
#markers=['o','s','D','^','<','v']
#assert len(colors) >= len(X)
xmin = ymin = sys.maxint
ymax=xmax=-sys.maxint
for key,col,mk,fl in zip(X,colors,markers,fillstyles):
if np.max(X[key]) > xmax:
xmax = np.max(X[key])
if np.min(rho_dict[key]) < ymin:
ymin = np.min(rho_dict[key])
if np.max(rho_dict[key]) > ymax:
ymax = np.max(rho_dict[key])
ax.errorbar(X[key],
rho_dict[key],
yerr=Err[key],
label='%2.1f-%2.1f'%(m_range[key][0],m_range[key][1]),
fmt='-o', markersize=8,color=col,marker=mk,markerfacecolor=fl,markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1
)
(ylo,yhi)=(10**np.floor(np.log10(ymin)),\
10**np.ceil(np.log10(ymax)))
xhi=10**np.ceil(np.log10(xmax))
ax.axis([xhi/(yhi/ylo),xhi,ylo,yhi])
#--- add major xticks
xmin=np.ceil(np.log10(ax.axis()[0]))
xmax=np.floor(np.log10(ax.axis()[1]))
nbin = xmax - xmin
ax.set_xticks(np.logspace(xmin,xmax,int(nbin)+1))
#--- put minor bins
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
#fig.legend(title='m',prop=fontP,loc='lower left',frameon=False)#, bbox_to_anchor=(0.4, .6))
#ax.tick_params(axis='y',left=True, right=True,which='both')
#ax.tick_params(axis='x',bottom=True, top=True,which='both')
DrawFrame(ax, (0.22,0.06),(0.14,0.06),0.01,LOG_X=True,LOG_Y=True)
#plt.tight_layout()
plt.savefig('%s/rho_r.png'%DIR_OUTPT_figs,dpi=300,bbox_inches='tight',transparent=True)
plt.show()
In [62]:
#--- choose bins for rho and lambda
#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi.reset_index(inplace=True,drop=True)
#--- remove pairs with r > cs.t
#cs = 3 * (24*3600) #--- km per day
df_trig = nij_trig.copy()
#--- parents' magnitudes
m_list = GetParentsMagnitudes( swarm_lohi, df_trig )
#--- number of daughters for each mother
df_nn = df_trig.groupby(by='Parent_id').count()['Event_id']
#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
#--- scattered plot
fig = plt.figure(figsize=(4,4),frameon=True)
ax = fig.add_subplot(111)
#ax.set_xlim(mc,swarm_lohi['magnitude'].max()+0.5)
m_min = min(m_list)-1.0e-6
m_max = max(m_list)+1.0e-6 #--- exclude the largest event
ax.set_xlim(m_min,np.ceil(m_max))
#ax.set_ylim(0.5,1000)
ax.plot( swarm_lohi['magnitude'].loc[df_nn.index],
df_nn,
'o',color='black',alpha=.1)
ax.plot( swarm_lohi['magnitude'].loc[df_zeros.index],
df_zeros,
'x',color='gray')
plt.yscale('log')
#plt.xlabel('$m$')
#plt.ylabel(r'$N_{as}$')
#--- choose m0<m<m1
dm = 0.7
m_min = mc-1.0e-6
m_max = 5.0 #swarm_lohi['magnitude'].max()+1.0e-6 #--- exclud the largest one
n = int( ( m_max - m_min ) / dm )
m = list(np.linspace(m_min,m_max,n))
#df_parents = df_trig.groupby(by='Parent_id').groups #--- parent events
#ds = swarm_lohi['magnitude'].loc[df_parents.keys()] #--- mag. of parent events
#m=AdjustBins(ds,m)
#m.pop(-2)
m+=[6.0]
Mlist=m
m_mean=np.array(np.zeros(len(m) - 1))
m_err=np.array(np.zeros(len(m) - 1))
N_as = np.array(np.zeros(len(m) - 1))
err_N = np.array(np.zeros(len(m) - 1))
for i in xrange( len(m) - 1 ):
try:
m0 = m[ i ]
m1 = m[ i + 1 ]
df, m_mean[ i ], m_err[i] = GetPairsWithSpecifiedParentMag((m0,m1),
swarm_lohi,df_trig) #--- get parent with m0<m<m1
df.reset_index( inplace = True, drop = True )
# pdb.set_trace()
#--- number of daughters for each mother
df_n = df.groupby(by='Parent_id').count()['Event_id']
#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
subset = swarm_lohi.loc[list_of_zero_afs]
list_of_zero_afs = subset[(subset['magnitude'] >= m0) & (subset['magnitude'] <= m1)].index
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
#--- mean
# print [m0,m1],list_of_zero_afs
# pdb.set_trace()
N_as[ i ] = 1.0*(df_n.sum()+ df_zeros.sum())/(len(df_n)+len(df_zeros))
err_N[ i ] = (N_as[ i ]/(len(df_n)+len(df_zeros)))**0.5 #--- poisson
# err_N[ i ]=np.std(df_n.tolist()+df_zeros.tolist())/(len(df_n)+len(df_zeros))**0.5
except:
# traceback.print_exc()
continue
#ax.plot(m, np.append(N_as,N_as[-1]), linestyle='-', drawstyle='steps-post',color='red')
ax.errorbar(m_mean,N_as,
yerr=err_N,xerr=m_err,
markerfacecolor='white',
fmt='o', markersize=8,color='red',markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1)
#--- add major xticks
xmin=np.ceil(ax.axis()[0])
xmax=np.floor(ax.axis()[1])
nbin = xmax - xmin
ax.set_xticks(np.linspace(xmin,xmax,int(nbin)+1))
DrawFrame(ax, (0.15,0.08),(0.14,0.06),0.1,LOG_Y=True)
try:
from scipy.optimize import curve_fit
def func(x,c,alpha):
return c*10**(alpha*x)
popt, pconv = curve_fit(func,m_mean[m_mean>=mc],N_as[m_mean>=mc])
mm = swarm_lohi['magnitude'].loc[df_zeros.index].tolist() + \
swarm_lohi['magnitude'].loc[df_nn.index].tolist() #swarm_lohi[swarm_lohi['magnitude']>=mc]['magnitude']
mm = np.array(mm)
mm=list(set(mm))
mm.sort()
mm=np.array(mm)
ax.plot(mm,func(mm,*popt),'-.g')
print popt
print pconv**0.5
except:
pass
#plt.savefig('%s/prod.png'%DIR_OUTPT_figs,dpi=150,bbox_inches='tight')
fig = plt.figure(figsize=(4,4),frameon=True)
ax = fig.add_subplot(111)
ax.plot(m, np.append(N_as,N_as[-1]), linestyle='-', drawstyle='steps-post',color='red')
ax.errorbar(m_mean,N_as,
yerr=err_N,xerr=m_err,
markerfacecolor='white',
fmt='o', markersize=8,color='red',markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1)
ax.set_xlim(m_min,np.ceil(m_max))
#ax.set_ylim(0.5,1000)
ax.set_yscale('log')
ax.set_xticks(np.linspace(xmin,xmax,int(nbin)+1))
plt.show()
In [63]:
#--- omori relation
import numpy as np
import math
from scipy.optimize import curve_fit
#def func(x,c,p,A):#,B):
# return A/(1+x/c)**p
def func(x,c,p,A,c2,p2):#,B):
return A/(1+(x/c)**p+(x/c2)**p2)
#def func(x,c,p):#,B):
# return (p-1)/c/(1+x/c)**p
#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi.reset_index(inplace=True,drop=True)
popt={}
pcov={}
perr = {}
dm = 0.7
m_min = mc-1.0e-6
m_max = 5.0 #swarm_lohi['magnitude'].max()+1.0e-6 #--- exclud the largest one
n = int( ( m_max - m_min ) / dm )
m = list(np.linspace(m_min,m_max,n))
#df_parents = df_trig.groupby(by='Parent_id').groups #--- parent events
#ds = swarm_lohi['magnitude'].loc[df_parents.keys()] #--- mag. of parent events
#m=AdjustBins(ds,m)
#m.pop(-2)
m+=[6.0]
Mlist=m
X={}
rho_dict={}
#m_range={}
Err={}
m_mean=np.array(np.zeros(len(m)-1))
Xedges={}
nrows=len(m)-1;ncols=2
fig=plt.figure(figsize=(ncols*4,nrows*4))
for i in xrange( len(m) - 1 ):
# if i == (len(m) - 2):
# pdb.set_trace()
m0 = m[ i ]
m1 = m[ i + 1 ]
df_trig,m_mean[ i ],junk = GetPairsWithSpecifiedParentMag((m0,m1),swarm_lohi,nij_trig) #--- get parent with m0<m<m1
df_trig.reset_index( inplace = True, drop = True )
print "# of mother events:%s"%len(df_trig.groupby(by="Parent_id"))
#--- number of daughters for each mother
df_n = df_trig.groupby(by='Parent_id').count()['Event_id']
#--- include daughters with zero children
list_of_zero_afs = set(df_trig['Event_id'])-set(df_trig['Parent_id'])
subset = swarm_lohi.loc[list_of_zero_afs]
list_of_zero_afs = subset[(subset['magnitude'] >= m0) & (subset['magnitude'] <= m1)].index
df_zeros = pd.Series(np.zeros(len(list_of_zero_afs)),index=list_of_zero_afs)
#--- mean
N_as = 1.0*(df_n.sum()+ df_zeros.sum())/(len(df_n)+len(df_zeros))
print 'N_as=%s'%N_as
#--- add distance & time column
# df_trig['R'] = rmat[df_trig['Parent_id'],df_trig['Event_id']]
# df_trig['T'] = tmat[df_trig['Parent_id'],df_trig['Event_id']]
# assert len ( df_trig[ df_trig['R'] == 0.0 ] ) == 0, '%s'%display( df_trig[ df_trig['R'] == 0.0 ] )
#--- remove pairs with r > cs.t
# cs = 3 * (24*3600) #--- km per day
# df_trig = RemovePair( df_trig, cs )
#--- rho plots
n_decades = int( math.ceil( log(df_trig['T'].max()/df_trig['T'].min(),10) ) )
nbin_per_decade = 4
nbins = nbin_per_decade * n_decades
bins = np.logspace(log(df_trig['T'].min(),10)-1.e-10,
log(df_trig['T'].max(),10)+1.e-10,nbins)
bins = AdjustBins(df_trig['T'],bins)
rho, xedges=np.histogram(df_trig['T'], #--- density
density=True,
bins=bins)
rho *= N_as #--- int rho = N_as
hist, xedges=np.histogram(df_trig['T'], #--- histogram
bins=bins)
x_bins, junk=np.histogram(df_trig['T'], #--- histogram
bins=bins,
weights = df_trig['T'])
#--- rm empty bins
x_bins = x_bins[ hist != 0 ]
rho = rho[ hist != 0 ]
hist = hist[ hist != 0 ]
#--- poisson err
err = rho/np.sqrt(hist)
x_bins /= hist
#---- fit
# fit = plfit.plfit( np.array(df_trig['T']) ) #--- clauset
# cval = fit[1] #--- completeness mag.
# pval = fit[0]
# [pval, gof] = plpva.plpva(np.array(df_trig['T']), cval, 'xmin', cval,'silent')
# print pval
#--- 2nd fit
try:
popt[i], pcov[i] = curve_fit(func, x_bins, rho,
p0=(1.0e-2,1.0,1e2,1e0,2.0),
bounds=((1.0e-6,0.1,1e-5,1e-6,0.1), (np.inf,7,np.inf,np.inf,7)), #c,p,A,c2,p2
sigma=err #rho**0.5 #err
# absolute_sigma = True,
# method='trf'#'dogbox'
)
perr[i] = np.sqrt(np.diag(pcov[i]))
except:
pass
#--- scattered plot
ax=fig.add_subplot(nrows,ncols,2*i+1)
ax.set_title(r'$%2.1f<m<%2.1f$'%(m0,m1))
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('$t$')
ax.set_ylabel('$r$')
#plt.xlim(1e-3,100)
#plt.ylim(1e-3,1e2)
ax.scatter(df_trig['T'],df_trig['R'],alpha=0.1)
#--- discretization effects
dt_min=df_trig['T'].min()
dt_max=df_trig['T'].max()
dr_min=df_trig['R'].min()
dr_max=df_trig['R'].max()
ax.plot([dt_min,dt_min],[dr_min,dr_max],
'-',color='blue')
ax.plot([dt_max,dt_max],[dr_min,dr_max],
'-',color='blue')
ax.plot([dt_min,dt_max],[dr_min,dr_min],
'-',color='blue')
ax.plot([dt_min,dt_max],[dr_max,dr_max],
'-',color='blue')
ax.set_ylim(dr_min/2,dr_max*2)
ax.plot([dt_min,dt_max],[dt_min*cs,dt_max*cs],
'-',color='red')
#--- rho
ax=fig.add_subplot(nrows,ncols,2*i+2)
ax.set_xlim(1e-4,1e2)
# ax.set_ylim(1e-3,1e3)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('$t$')
ax.set_ylabel('$\lambda(t)$')
ax.errorbar(x_bins,rho,
yerr=err,
linestyle='None',
fmt='-o', markersize=8
)
#--- plot the fit
# plt.plot(x_bins,
# (pval-1)/cval/(1+x_bins/cval)**pval,
# label='fit: p=%2.1f, c=%4.3f' % (pval,cval)
# )
#--- 2nd fit
try:
ax.plot(x_bins, func(x_bins, *popt[i]),'r-',label='fit: p0=%3.2f, p1=%3.2f' % (popt[i][1],popt[i][4]))
except:
pass
fig.tight_layout()
# ax.legend(prop=fontP)#bbox_to_anchor=(1.5, 0.5))
X[i] = x_bins[:]
Xedges[i]=xedges
rho_dict[i]=rho[:]
# m_range[i]=np.array([m0,m1])
Err[i]=err[:]
#fig.savefig('%s/rho_t.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight')
plt.show()
In [64]:
import traceback
try:
KEYS=popt.keys()
KEYS.sort()
M_Mean=[]
pval=[]
cval = []
cerr = []
aval = []
aerr = []
mm=[]
pperr=[]
fig = plt.figure(figsize=(12,4))
ax=fig.add_subplot(131)
for ikey in KEYS: #--- exclude the last point
pval.append( min(popt[ikey][1],popt[ikey][4]) )
cval.append( min(popt[ikey][0],popt[ikey][3]) )
aval.append( popt[ikey][2] )
aerr.append(perr[ikey][2])
index = min([popt[ikey][1],1],[popt[ikey][4],4])[1]
pperr.append(perr[ikey][index])
cerr.append(perr[ikey][index-1])
M_Mean.append(m_mean[ikey])
ax.errorbar(M_Mean,pval,yerr=pperr,
marker='o',color='black')
#ax.set_xlim(0,4)
ax.set_ylim(0,3)
ax.set_xlabel('$m$')
ax.set_ylabel('$p$')
ax.set_title('nth=%2.1e'%n_thresh)
#ax.plot([0,3],[1.29,1.29],'--r')
#ax.plot([0,3],[1.37,1.37],'--r')
#--- plot c
def func(x,a,b):
# return a*10**(x*b)
return 10**(b+x*a)
# poptt, pcovv = curve_fit(func, M_Mean,cval, sigma=cerr)
poptt = np.polyfit(M_Mean,np.log10(cval),1)
alpha_c = poptt[0]
# poptt[1] = 10**poptt[0]
# print poptt
ax2=fig.add_subplot(132)
# ax2.set_ylim(1e-4,1e-1)
#ax2.set_xlim(0,4)
ax2.set_yscale('log')
ax2.set_ylabel('$c$')
ax2.set_xlabel('$m$')
ax2.errorbar(M_Mean,cval,yerr=cerr,
marker='o',color='black')
ax2.plot(M_Mean,func(np.array(M_Mean),*poptt),'--r')
# #--- plot A
# poptt, pcovv = curve_fit(func, M_Mean,aval, sigma=aerr)
# ax3=fig.add_subplot(133)
# #ax3.set_ylim(1e2,1e4)
# #ax3.set_xlim(0,4)
# ax3.set_yscale('log')
# ax3.set_ylabel('$A$')
# ax3.set_xlabel('$m$')
# ax3.errorbar(M_Mean,aval,yerr=aerr,
# marker='o',color='black')
# ax3.plot(m_mean,func(M_Mean,*poptt),'--r')
# print poptt
fig.tight_layout()
#plt.savefig('%s/p.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight')
except:
traceback.print_exc()
pass
plt.show()
In [65]:
#--- omori relation
plt.figure(figsize=(16.5,5))
#--- unrescaled
plt.subplot(131)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('$t$')
plt.ylabel('$\lambda(t)$')
#plt.xlim(1e-6,1e1)
#plt.ylim(1e0,1e7)
KEYS = X.keys()
KEYS.sort()
for key, index in zip(KEYS,xrange(100)):
plt.errorbar(X[key],rho_dict[key],yerr=Err[key],
# label='%2.1f-%2.1f'%(m_range[key][0],m_range[key][1]),
fmt='-o', markersize=8
)
#plt.legend(title='m',fontsize='small',loc='upper right',bbox_to_anchor=(0.4, .65),frameon=False)
xdata = []
ydata = []
XXdata = []
YYdata = []
y_err = []
#--- rescaled
#alpha_c = 0.0 #0.5
#alpha_p = 0.43 #--- productivity
alpha_y = alpha_c - alpha_p
plt.subplot(132)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'$t/10^{%3.2f\times m}$'%alpha_c)
plt.ylabel(r'$\lambda(t) \times 10^{(%3.2f-%3.2f)m}$'%(alpha_c,alpha_p))
#plt.ylabel(r'$\lambda(t) \times t^p$')
#plt.xlim(1e-6,1e0)
#plt.ylim(1e-3,1e5)
for key, index in zip(KEYS,xrange(100)):
if key > 3:
continue
xdata += list(X[key]/10**(alpha_c*m_mean[key]))
ydata += list(rho_dict[key]*10**(alpha_y*m_mean[key]))
y_err += list(Err[key]*10**(alpha_y*m_mean[key]))
XXdata += list(X[key])
YYdata += list(rho_dict[key]*10**(-alpha_p*m_mean[key]))
plt.errorbar(X[key]/10**(alpha_c*m_mean[key]),
rho_dict[key]*10**(alpha_y*m_mean[key]),
yerr=Err[key]*10**(alpha_y*m_mean[key]),
# label='%2.1f-%2.1f'%(m_range[key][0],m_range[key][1]),
label='%2.1f'%(m_mean[key]),
fmt='-o', markersize=8
)
plt.legend(title='m',fontsize='small',loc='upper right',bbox_to_anchor=(0.4, .65),frameon=False)
#--- rescaled by t^-p
p = 1
plt.subplot(133)
plt.xscale('log')
plt.yscale('log')
plt.xlabel(r'$t/10^{%3.2f\times m}$'%alpha_c)
plt.ylabel(r'$\lambda(t) \times 10^{(%3.2f-%3.2f)m} \times t^{%3.2f}$'%(alpha_c,alpha_p,p))
#plt.ylabel(r'$\lambda(t) \times t^p$')
#plt.xlim(1e-6,1e0)
#plt.ylim(1e-5,1e-3)
for key, index in zip(KEYS,xrange(100)):
x_new = X[key]/10**(alpha_c*m_mean[key])
y_new = rho_dict[key]*10**(alpha_y*m_mean[key])
y2_new = Err[key]*10**(alpha_y*m_mean[key])
plt.errorbar(x_new,
y_new*x_new**p,
yerr=y2_new*x_new**p,
fmt='-o', markersize=8
)
plt.tight_layout()
###plt.savefig('%s/rho_t_mult.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight')
plt.show()
In [66]:
#--- fitting
from scipy.optimize import curve_fit
def func(x,c,p,A):#,c2,p2):#,B):
return A/(1+(x/c)**p) #+(x/c2)**p2)
#def func(x,b,a):
# return a*x**b
#tc=1e-3
xdata = np.array(xdata)
xxdata = xdata[:]#[xdata>tc]
ydata = np.array(ydata)
yydata = ydata[:] #xdata>tc]
Y_err=y_err
#
ss=zip(xxdata,yydata,Y_err)
ss.sort()
xxdata=[i[0]for i in ss]
yydata=[i[1]for i in ss]
Y_err = [i[2] for i in ss]
popt, pcov = curve_fit(func, xxdata, yydata,
p0=(1.0e-2,1.0,1e2),#,1e0,2.0),
bounds=((1.0e-6,0.1,1e-5), (np.inf,7,np.inf)), #c,p,A,c2,p2
sigma=Y_err)
p = -popt[1]
rho_pref = func(np.min(xxdata),*popt)
tlo = popt[0]
#--- linear
#popt = np.polyfit(np.log10(xxdata),np.log10(yydata), 1)
#popt[1] = 10**(popt[1])
#p = popt[0]
#rho_pref = popt[1]
#--- plot
fig = plt.figure(figsize=(4,4))
ax=fig.add_subplot(111)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('$t$')
ax.set_ylabel(r'$\lambda(t)$')
#ax.set_xlim(2,6)
#ax.set_ylim(1e-4,1e2)
ax.scatter(xdata,ydata)
ax.scatter(xdata,1e2*ydata*xdata**(-p))
ax.scatter([popt[0]],func([popt[0]],*popt),
marker='x',c='black',s=160)
#--- fit
ax.plot(xxdata,func(xxdata,*popt),
'-.r',zorder=3,
# label=r'$\lambda(t)\propto t^{%3.2f}$'%popt[0])
label=r'$\lambda(t)\propto t^{%3.2f}$'%(p))
ax.legend(fontsize='xx-small',frameon=False)
#ax.set_title('q=%d%%,df=%s'%(quantile*100,Df),fontsize='xx-small')
plt.savefig('%s/omori_1st.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight',transparent=True)
#--- save p
np.save('omori_exp1st',p)
plt.show()
In [67]:
alpha_arr=[]
for irun in xrange(16):
path = 'runYuha/Run%s'%irun
a=np.load('%s/omori_exp1st.npy'%path)
alpha_arr.append(a.item())
#--- load image
display(Image(filename='%s/omori_1st.png'%path))
df=pd.DataFrame({'mean':[np.mean(alpha_arr)], 'std_dev':np.std(alpha_arr),
'std_err':[np.std(alpha_arr)/len(alpha_arr)**.5]},
index=[r'$\sigma$'])
display(df)
df['mean']+2*df['std_dev'],df['mean']-2*df['std_dev']
In [68]:
#--- fitting
from scipy.optimize import curve_fit
def func(x,b,a):
return a*x**b
#--- 2nd regime
xdata = np.array(xdata)
xxdata = xdata[xdata>tlo]
#
ydata = np.array(ydata)
yydata = ydata[xdata>tlo]
#
y_err = np.array(y_err)
Y_err=y_err[xdata>tlo]
#
ss=zip(xxdata,yydata,Y_err)
ss.sort()
xxdata=[i[0]for i in ss]
yydata=[i[1]for i in ss]
Y_err = [i[2] for i in ss]
#--- linear
popt = np.polyfit(np.log10(xxdata),np.log10(yydata), 1,
w=np.array(yydata)/np.array(Y_err))
popt[1] = 10**(popt[1])
p = popt[0]
rho_pref = popt[1]
#--- plot
fig = plt.figure(figsize=(4,4))
ax=fig.add_subplot(111)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('$t$')
ax.set_ylabel(r'$\lambda(t)$')
#ax.set_xlim(2,6)
#ax.set_ylim(1e-4,1e2)
ax.scatter(xdata,ydata)
ax.scatter(xdata,1e1*ydata*xdata**(-p))
#ax.scatter([popt[0]],func(np.array([popt[0]]),*popt), #--- size
# marker='x',c='black',s=100)
#--- fit
ax.plot(xxdata,func(xxdata,*popt),
'-.r',zorder=3,
label=r'$\lambda(t)\propto t^{%3.2f}$'%(p))
ax.legend(fontsize='xx-small',frameon=False)
#ax.set_title('q=%d%%,df=%s'%(quantile*100,Df),fontsize='xx-small')
plt.savefig('%s/omori_2nd.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight',transparent=True)
#--- save p
np.save('omori_exp2nd',p)
plt.show()
In [69]:
alpha_arr=[]
for irun in xrange(16):
path = 'runYuha/Run%s'%irun
a=np.load('%s/omori_exp2nd.npy'%path)
alpha_arr.append(a.item())
#--- load image
display(Image(filename='%s/omori_2nd.png'%path))
df=pd.DataFrame({'mean':[np.mean(alpha_arr)], 'std_dev':np.std(alpha_arr),
'std_err':[np.std(alpha_arr)/len(alpha_arr)**.5]},
index=[r'$p$'])
display(df)
df['mean']+2*df['std_dev'],df['mean']-2*df['std_dev']
In [70]:
# #--- fitting
# from scipy.optimize import curve_fit
# def func(x,c,p,A,c2,p2):#,B):
# return A/(1+(x/c)**p+(x/c2)**p2)
# #def func(x,b,a):
# # return a*x**b
# #tc=1e-3
# xdata = np.array(xdata)
# xxdata = xdata[:]#[xdata>tc]
# ydata = np.array(ydata)
# yydata = ydata[:] #xdata>tc]
# Y_err=y_err
# #
# ss=zip(xxdata,yydata,Y_err)
# ss.sort()
# xxdata=[i[0]for i in ss]
# yydata=[i[1]for i in ss]
# Y_err = [i[2] for i in ss]
# popt, pcov = curve_fit(func, xxdata, yydata,
# p0=(1.0e-2,1.0,1e2,1e0,2.0),
# bounds=((1.0e-6,0.1,1e-5,1e-6,0.1), (np.inf,7,np.inf,np.inf,7)), #c,p,A,c2,p2
# sigma=Y_err)
# p = -popt[4]
# rho_pref = func(np.min(xxdata),*popt)
# #--- linear
# #popt = np.polyfit(np.log10(xxdata),np.log10(yydata), 1)
# #popt[1] = 10**(popt[1])
# #p = popt[0]
# #rho_pref = popt[1]
# #--- plot
# fig = plt.figure(figsize=(4,4))
# ax=fig.add_subplot(111)
# ax.set_xscale('log')
# ax.set_yscale('log')
# ax.set_xlabel('$t$')
# ax.set_ylabel(r'$\lambda(t)$')
# #ax.set_xlim(2,6)
# #ax.set_ylim(1e-4,1e2)
# ax.scatter(xdata,ydata)
# #--- fit
# ax.plot(xxdata,func(xxdata,*popt),
# '-.r',zorder=3,
# # label=r'$\lambda(t)\propto t^{%3.2f}$'%popt[0])
# label=r'$\lambda(t)\propto t^{%3.2f,%3.2f}$'%(-popt[4],-popt[1]))
# ax.legend(fontsize='xx-small',frameon=False)
# #ax.set_title('q=%d%%,df=%s'%(quantile*100,Df),fontsize='xx-small')
# ###plt.savefig('%s/rupture_mag.png'%DIR_OUTPT_figs,dpi=75,bbox_inches='tight',transparent=True)
# plt.show()
In [71]:
#--- omori (paper version)
#--- rescaled
fig, ax = plt.subplots(figsize=(4,4))
#ax.axis([1/(20/3e-3),1,3e-3,20])
#ax.loglog()
#colors = ['black','red','green','blue','cyan','brown']
#fillstyles=['white',None,'white',None,'white',None,'white']
#markers=['o','s','D','^','<','v']
#assert len(colors) >= len(X)
for key,col,mk,fl in zip(X,colors,markers,fillstyles):
ax.errorbar(X[key]/10**(alpha_c*m_mean[key]),
rho_dict[key]*10**(alpha_y*m_mean[key]),
yerr=Err[key]*10**(alpha_y*m_mean[key]),
# label='%2.1f-%2.1f'%(m_range[key][0],m_range[key][1]),
fmt='-o', markersize=8,color=col,marker=mk,markerfacecolor=fl,markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1
)
(ylo,yhi)=(10**np.floor(np.log10(np.min(ydata))),\
10**np.ceil(np.log10(np.max(ydata))))
xhi=10**np.ceil(np.log10(np.max(xdata)))
ax.axis([xhi/(yhi/ylo),xhi,ylo,yhi])
ax.loglog()
##--- add major xticks
xmin=np.ceil(np.log10(ax.axis()[0]))
xmax=np.floor(np.log10(ax.axis()[1]))
nbin = xmax - xmin
ax.set_xticks(np.logspace(xmin,xmax,int(nbin)+1))
#--- add major yticks
ymin=np.ceil(np.log10(ax.axis()[2]))
ymax=np.floor(np.log10(ax.axis()[3]))
nbin = ymax - ymin
ax.set_yticks(np.logspace(ymin,ymax,int(nbin)+1))
#--- put minor bins
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
#--- put minor bins y
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
#fig.legend(title='m',prop=fontP,loc='upper right',frameon=False)#, bbox_to_anchor=(0.4, .6))
ax.tick_params(axis='y',left=True, right=True,which='both')
ax.tick_params(axis='x',bottom=True, top=True,which='both')
ax.set_xlabel('$t(day)$')
ax.set_ylabel('$\lambda(t)/10^{%3.2fm}$'%alpha_p)
#--- plot power-law
(xlo,xhi)=(10**xmin,10**(xmin+6)) #(0.8*15e-3,0.8*.6)
yhi=func(xlo,*popt)*10 #rho_pref #*xmin**p/10
ax.plot([xlo,xhi],[yhi, yhi* (xhi/xlo)**p],'r-.',
label=r'$\lambda(t)\propto t^{%3.2f}$'%(p),
linewidth=1)
ax.legend(fontsize='x-small',frameon=False)
ax.set_xlabel(r'$t$')
#ax.set_ylabel(r'$\rho(r)\times r_{rup}$')
ax.set_title('mc=%s,~b=%s'%(mc,bval),fontsize='xx-small')
DrawFrame(ax, (0.22,0.06),(0.14,0.06),0.01,LOG_X=True,LOG_Y=True)
#plt.tight_layout()
plt.savefig('%s/rho_t2nd.png'%DIR_OUTPT_figs,dpi=75*2,bbox_inches='tight')
plt.show()
In [72]:
#--- omori (paper version)
import re
# create a definition for the short hyphen
matplotlib.rcParams["text.latex.preamble"].append(r'\mathchardef\mhyphen="2D')
class MyLogFormatter(matplotlib.ticker.LogFormatterMathtext):
def __call__(self, x, pos=None):
# call the original LogFormatter
rv = matplotlib.ticker.LogFormatterMathtext.__call__(self, x, pos)
# check if we really use TeX
if matplotlib.rcParams["text.usetex"]:
# if we have the string ^{- there is a negative exponent
# where the minus sign is replaced by the short hyphen
rv = re.sub(r'\^\{-', r'^{\mhyphen', rv)
return rv
#--- rescaled
fig, ax = plt.subplots(figsize=(4,4))
#ax.axis([1/(20/3e-3),1,3e-3,20])
#ax.loglog()
#colors = ['black','red','green','blue','cyan','brown']
#fillstyles=['white',None,'white',None,'white',None,'white']
#markers=['o','s','D','^','<','v']
#assert len(colors) >= len(X)
for key,col,mk,fl in zip(X,colors,markers,fillstyles):
ax.errorbar(X[key], #/10**(alpha_c*m_mean[key]),
rho_dict[key]*10**(-alpha_p*m_mean[key]),
yerr=Err[key]*10**(-alpha_p*m_mean[key]),
# label='%2.1f-%2.1f'%(m_range[key][0],m_range[key][1]),
fmt='-o', markersize=8,color=col,marker=mk,markerfacecolor=fl,markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1
)
(ylo,yhi)=(10**(np.floor(np.log10(np.min(YYdata)))-0),\
# 10**np.ceil(np.log10(np.max(YYdata))))
np.max(YYdata))
#(xlo,xhi)=(10**np.floor(np.log10(np.min(xdata))),\
# 10**np.ceil(np.log10(np.max(xdata))))
xc=(np.min(XXdata)*np.max(XXdata))**.5
xlo = xc/(yhi/ylo)**.5
xhi = xc*(yhi/ylo)**.5
ax.axis([xlo,xhi,ylo,yhi])
ax.loglog()
##--- add major xticks
xmin=np.ceil(np.log10(ax.axis()[0]))
xmax=np.floor(np.log10(ax.axis()[1]))
nbin = xmax - xmin
ax.set_xticks(np.logspace(xmin,xmax,int(nbin)+1))
#--- add major yticks
ymin=np.ceil(np.log10(ax.axis()[2]))
ymax=np.floor(np.log10(ax.axis()[3]))
nbin = ymax - ymin
ax.set_yticks(np.logspace(ymin,ymax,int(nbin)+1))
#--- put minor bins
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
ax.xaxis.set_major_formatter(MyLogFormatter()) #--- minus sign too long
#--- put minor bins y
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
ax.xaxis.set_major_formatter(MyLogFormatter()) #--- minus sign too long
#fig.legend(title='m',prop=fontP,loc='upper right',frameon=False)#, bbox_to_anchor=(0.4, .6))
ax.tick_params(axis='y',left=True, right=True,which='both')
ax.tick_params(axis='x',bottom=True, top=True,which='both')
#ax.set_xlabel('$t(day)$')
#ax.set_ylabel('$\lambda(t)/10^{%3.2fm}$'%alpha_p)
#--- plot power-law
tc=tlo
(xlo,xhi)=(tc,1e7*tc) #(0.8*15e-3,0.8*.6)
yhi=rho_pref*tc**p*10*10
ax.plot([xlo,xhi],[yhi, yhi* (xhi/xlo)**p],'r-.',
# label=r'$\lambda(t)\propto t^{%3.2f}$'%(p),
linewidth=1)
#ax.legend(fontsize='x-small',frameon=False)
#ax.set_xlabel(r'$t$')
#ax.set_ylabel(r'$\rho(r)\times r_{rup}$')
#ax.set_title('q=%d%%,df=%s'%(quantile*100,Df),fontsize='xx-small')
DrawFrame(ax, (0.22,0.06),(0.14,0.06),0.1,LOG_X=True,LOG_Y=True)
#plt.tight_layout()
plt.savefig('%s/rho_t2nd.png'%DIR_OUTPT_figs,dpi=300,bbox_inches='tight')
plt.show()
In [73]:
#--- density plots: rho(r) (paper version 2nd)
#plt.style.use('seaborn')
#plt.rcParams.update({'lines.markeredgewidth': 1})
def returnCMtoInch(x,y):
return (x/2.54,y/2.54)
fig, ax = plt.subplots(figsize=(4,4))
#---
ax.spines['right'].set_visible(False) #--- half open
ax.spines['top'].set_visible(False)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
#---
#ax.axis([1e-3,1e2,1e-4,1e1])
ax.loglog()
#ax.set_xticks(np.logspace(-1,3,4))
#plt.xlabel(r'$R=r/10^{%s m}$'%alpha)
#plt.ylabel(r'$P(R)=P(r)\times 10^{%s m}$'%alpha)
#plt.xlim(1e-5,1e0)
#plt.ylim(1e-5,1e0)
#colors = ['black','red','green','blue','cyan','brown']
#fillstyles=['white',None,'white',None,'white',None,'white']
#markers=['o','s','D','^','<','v']
#assert len(colors) >= len(X)
xmin = ymin = sys.maxint
ymax=xmax=-sys.maxint
for key,col,mk,fl in zip(X,colors,markers,fillstyles):
if np.max(X[key]) > xmax:
xmax = np.max(X[key])
if np.min(X[key]) < xmin:
xmin = np.min(X[key])
if np.min(rho_dict[key]) < ymin:
ymin = np.min(rho_dict[key])
if np.max(rho_dict[key]) > ymax:
ymax = np.max(rho_dict[key])
ax.errorbar(X[key],
rho_dict[key],
yerr=Err[key],
label='%2.1f-%2.1f'%(m_range[key][0],m_range[key][1]),
fmt='-o', markersize=8,color=col,marker=mk,markerfacecolor=fl,markeredgewidth=0.7,
linewidth=.5,
barsabove=None,capsize=5,capthick=1,elinewidth=1
)
(ylo,yhi)=(10**np.floor(np.log10(ymin)),\
# 10**np.ceil(np.log10(ymax)))
ymax)
xc=(xmin*xmax)**.5
xlo = xc/(yhi/ylo)**.5
xhi = xc*(yhi/ylo)**.5
ax.axis([xlo,xhi,ylo,yhi])
#--- add major xticks
xmin=np.ceil(np.log10(ax.axis()[0]))
xmax=np.floor(np.log10(ax.axis()[1]))
nbin = xmax - xmin
ax.set_xticks(np.logspace(xmin,xmax,int(nbin)+1))
#--- add major yticks
ymin=np.ceil(np.log10(ax.axis()[2]))
ymax=np.floor(np.log10(ax.axis()[3]))
nbin = ymax - ymin
ax.set_yticks(np.logspace(ymin,ymax,int(nbin)+1))
#--- put minor bins
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.xaxis.set_minor_locator(locmin)
ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
ax.xaxis.set_major_formatter(MyLogFormatter()) #--- minus sign too long
#--- put minor bins y
locmin = matplotlib.ticker.LogLocator(base=10.0,subs=(0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1),numticks=12)
ax.yaxis.set_minor_locator(locmin)
ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())
ax.yaxis.set_major_formatter(MyLogFormatter())
#fig.legend(title='m',prop=fontP,loc='lower left',frameon=False)#, bbox_to_anchor=(0.4, .6))
#ax.tick_params(axis='y',left=True, right=True,which='both')
#ax.tick_params(axis='x',bottom=True, top=True,which='both')
DrawFrame(ax, (0.22,0.06),(0.14,0.06),0.01,LOG_X=True,LOG_Y=True)
#plt.tight_layout()
plt.savefig('%s/rho_t.png'%DIR_OUTPT_figs,dpi=300,bbox_inches='tight',transparent=True)
plt.show()
In [74]:
#--- construct trees
from math import *
class tree():
leafDepth = 0
n_leaves = 0
g_netx = None
leaf_generation = {} #--- key: leaf id; val: generation
nsize = 1 #--- cluster size
rtm_dic = {}
pairs = []
#---
def __init__( self, data, coord, coord_parent, n_generation ):
self.data = data #--- id
self.coord = coord
self.coord_parent = coord_parent
self.n_daughters = 0
self.daughter_obj = {}
self.n_generation = n_generation
#---
def insert( self, data, parent_id, coord, coord_parent, n_generation ):
if parent_id == self.data: #--- parent exists
self.daughter_obj[ self.n_daughters ] = tree( data, coord, coord_parent, n_generation )
self.n_daughters += 1
return
else:
for keys in self.daughter_obj:
self.daughter_obj[ keys ].insert( data, parent_id, coord, coord_parent, n_generation )
#---
def sprint( self, df ):
for keys in self.daughter_obj:
tree.g_netx.add_node(self.daughter_obj[ keys ].data)
tree.g_netx.add_edge( self.data, self.daughter_obj[ keys ].data)
self.daughter_obj[ keys ].sprint( df )
#---
def getLeaves( self ):
if len( self.daughter_obj ) == 0: #--- it's a leaf
tree.leaf_generation[ self.data ] = self.n_generation
# print tree.leafDepth, self.data, self.n_generation
else:
for keys in self.daughter_obj: #--- search in daughters
self.daughter_obj[ keys ].getLeaves()
#---
def getSize( self ):
for keys in self.daughter_obj: #--- search in daughters
tree.nsize += 1
# print 'id=%s,size=%s'%(self.daughter_obj[ keys ].data, tree.nsize)
self.daughter_obj[ keys ].getSize()
def get( self, title ):
pass
# tree.c.writePDFfile( title )
def get_coords( self ): #--- rtm_dic['node_id'] = [r,t,m]
for keys in self.daughter_obj:
child_object = self.daughter_obj[keys]
cords = child_object.coord
id_node = child_object.data
tree.rtm_dic[ id_node ] = cords
child_object.get_coords()
def plot(self,ax):
# print 'id=%s'%(self.data), 'generation=%s'%(self.n_generation)
for keys in self.daughter_obj:
child_object = self.daughter_obj[keys]
cords = child_object.coord
tree.pairs.append(child_object.data) #[self.coord, child_object.coord])
# ax.annotate("", tuple(cords[1][:2]), xytext=tuple(self.coord[1][:2]),
# textcoords='data', annotation_clip=False,
# arrowprops=dict(arrowstyle="-|>,head_width=.1,head_length=.4",color="b",linewidth="0.3"))
child_object.plot(ax)
def nodeExist(self, event_id):
if self.data == event_id:
# print 'yes'
return True
for keys in self.daughter_obj:
child_object = self.daughter_obj[keys]
if child_object.nodeExist( event_id ):
return True
return False
pref = 1.0 / (24.0*60*60)
#--------------------
#----- subset (the same as trig. analysis)
#--------------------
swarm_lohi = DataFrameSubSet( swarm,
'date',
( myTlohi.lo,myTlohi.hi ) )
swarm_lohi = DataFrameSubSet( swarm_lohi, #--- complete catalog
'magnitude',
( mc, sys.maxint ) )
swarm_lohi = GetCartesian( swarm_lohi )
swarm_lohi.reset_index(inplace=True,drop=True)
NoTriggerList = set(nij_trig['Parent_id'])-set(nij_trig['Event_id']) #--- list of roots
#---
#--- include parent and entire generation
#---
CLUSTER={}
clusterID = 0
stree = {}
magScale = 1.0
for rootID in NoTriggerList: #--- roots: events with no parents
tree_level = 0
mainShockList = [ rootID ]
CLUSTER[clusterID]=[mainShockList[0]]
t_root = swarm_lohi.loc[ rootID ]['date']
r_root = swarm_lohi.loc[ rootID ][['x(km)','y(km)','z(km)']].tolist()
m_root = swarm_lohi.loc[ rootID ]['magnitude']
stree[ clusterID ] = tree( rootID, [ t_root, r_root, m_root ],
[ t_root, r_root, m_root ],
tree_level )
while len(mainShockList): #--- from root down to leaves
newMainShockList = []
for mainShockID, counter in zip( mainShockList, xrange( sys.maxint ) ):
try:
newMainShockList += d_trig[ mainShockID ] #--- aftershocks become new
except:
# traceback.print_exc()
continue
for aftershockID in d_trig[ mainShockID ]:
CLUSTER[ clusterID ].append(aftershockID)
# dt = (swarm_lohi.loc[ aftershockID ]['date']-swarm_lohi.loc[ rootID ]['date']).total_seconds() * pref
t_afs = swarm_lohi.loc[ aftershockID ]['date']
# -swarm_lohi.loc[ rootID ]['date']).total_seconds() * pref
r_afs = swarm_lohi.loc[aftershockID][['x(km)','y(km)','z(km)']].tolist()
m_afs = swarm_lohi.loc[aftershockID]['magnitude']
t_parent = swarm_lohi.loc[ mainShockID ]['date']
r_parent = swarm_lohi.loc[mainShockID][['x(km)','y(km)','z(km)']].tolist()
m_parent = swarm_lohi.loc[mainShockID]['magnitude']
stree[ clusterID ].insert( aftershockID, mainShockID,
[t_afs, r_afs, m_afs],
[t_parent, r_parent, m_parent],
tree_level + 1 )
mainShockList = [ i for i in newMainShockList ]
tree_level += 1
assert len(mainShockList) == 0 #--- leaf
clusterID+=1
sl=[]
for keys in stree:
sl.append([stree[keys].n_daughters,stree[keys].data,keys])
sl.sort(reverse=True)
sl[:3]
Out[74]:
In [75]:
#--- structure of the trig. cascade
import networkx as nx
#import treeDataStruct as tds
from itertools import count
import random as rnd
for item, index in zip(sl[:6],xrange(sys.maxint)):
sid = item[2]
root_id = item[1]
print 'root_id=',root_id, 'm=', swarm_lohi['magnitude'].iloc[root_id]
tree.g_netx = None
tree.g_netx = nx.DiGraph(directed=True)
tree.g_netx.add_node(stree[ sid ].data,color='b')
stree[ sid ].sprint(swarm_lohi )
# read edgelist from grid.edgelist
nx.write_edgelist(tree.g_netx, path="grid.edgelist", delimiter=":")
H = nx.read_edgelist(path="grid.edgelist", delimiter=":") #,create_using=nx.DiGraph())
# nx.draw(H,with_labels=True,arrows = True)#pos=nx.spring_layout(H))
# nx.draw_networkx(tree.g_netx, arrows=True, with_labels=True,node_color='b') #, **kwds)
pos = nx.spring_layout(H)
nodes=H.nodes
colors = ['red' for n in nodes]
sizes=[3**swarm_lohi.loc[ int(n) ]['magnitude'] for n in nodes]
for n,indx in zip(nodes,xrange(sys.maxint)):
if int(n) == root_id:
colors[indx] = 'gray'
ec = nx.draw_networkx_edges(H, pos, alpha=0.2,arrows = True)
nc = nx.draw_networkx_nodes(H, pos,node_size=sizes, node_color=colors,cmap=plt.cm.seismic) #, nodelist=nodes, ,
# with_labels=False, node_size=100, cmap=plt.cm.jet)
plt.axis('off')
#plt.savefig('%s/topology.%s.png'%(DIR_OUTPT_figs,index),dpi=75)#,bbox_inches='tight')
plt.show()
size vs. depth¶
In [76]:
mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings
rcc('text', usetex=True)
font = {'size' : 26}
matplotlib.rc('font', **font)
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.set_xscale('linear')
ax.set_yscale('log')
#ax.set_xlabel('mean depth')
#ax.set_ylabel('size')
ax.axis([0,7.0,1,1e3])
ax.set_yscale('log')
#--- add major xticks
xmin=np.ceil(ax.axis()[0])
xmax=np.floor(ax.axis()[1])
nbin = xmax - xmin
ax.set_xticks(np.linspace(ax.axis()[0],ax.axis()[1],int(nbin)+1))
#--- add major yticks
#ymin=np.ceil(np.log10(ax.axis()[2]))
#ymax=np.floor(np.log10(ax.axis()[3]))
#nbin = ymax - ymin
#ax.set_yticks(np.logspace(ymin,ymax,int(nbin)+1))
MD=[]
S=[]
M=[]
root=[]
for item in sl:
sid = item[2]
root_id = item[1]
root.append(root_id)
tree.leaf_generation={} #--- initialize
# stree[ sid ].initialize()
stree[ sid ].getLeaves() #--- call function
d = stree[ sid ].leaf_generation
# print 'd[%s]=%s'%(sid,d)
mean_depth = np.mean( d.values() ) #np.mean( d.values() ) #1.0 * tds.tree.leafDepth / tds.tree.n_leaves
#--- size
tree.nsize = 1 #--- initialize
stree[ sid ].getSize() #--- call function
MD.append(mean_depth)
S.append(stree[ sid ].nsize)
M.append(swarm_lohi.loc[stree[ sid ].data]['magnitude'])
exp = 3.2
S=np.array(S)
MD = np.array(MD)
M=np.array(M)
root=np.array(root)
scatter=ax.scatter(MD, S,
marker='o',color='gray',edgecolors='none',alpha=0.4,
s=exp**np.array(M))
#--- color specific clusters
cond = np.all([S>=100,MD<2],axis=0)
scatter=ax.scatter(MD[cond], S[cond],
marker='o',color='red',edgecolors='none',#alpha=0.5,
s=exp**np.array(M[cond]))
#print root[10<=S], S[10<=S]
print 'roots=',root[cond]
print 'size',S[cond]
print 'depth',MD[cond]
# #--- legend
# M=np.array(M)
# m_min=M.min()
# m_max=M.max()
# m_mean=(m_min+m_max)/2 #M.mean()
# l1 = plt.scatter([-2],[-2], s=exp**m_min, c='gray',edgecolors='none')
# l2 = plt.scatter([-3],[-3], s=exp**m_mean,c='gray', edgecolors='none')
# l3 = plt.scatter([-4],[-4], s=exp**m_max,c='gray', edgecolors='none')
# labels = ["%2.1f"%m_min, "%2.1f"%m_mean, "%2.1f"%m_max]
# leg = plt.legend([l1, l2, l3], labels, ncol=3,fontsize=16, frameon=False,
# handlelength=0, borderpad = 0, labelspacing=0,columnspacing=2,prop=None,
# handletextpad=1, scatterpoints = 1,loc='upper right', bbox_to_anchor=(.92, 0.95))
#--- frame
DrawFrame(ax, (0.24,0.08),(0.14,0.06),0.1,LOG_Y=True)
#--- save
#plt.savefig('%s/topology.png'%DIR_OUTPT_figs,dpi=300,bbox_inches='tight')
mpl.rcParams.update(mpl.rcParamsDefault)
warnings.filterwarnings('ignore') #--- get rid of warnings
rcc('text', usetex=True)
font = {'size' : 26}
matplotlib.rc('font', **font)
In [77]:
x=(-97,-96.5)
y=(35.8,36)
swarm_lohi.sort_values(by='magnitude',ascending=False)[(swarm_lohi['longitude']>x[0]) & (swarm_lohi['longitude']<x[1]) & (swarm_lohi['latitude']>y[0]) & (swarm_lohi['latitude']<y[1])].head(n=10)
Out[77]:
In [78]:
#--- a large event belongs to what cluster
event_id = 13698 #91 #2917
FOUND = False
for item in sl:
cls_id = item[2]
root_id = item[1]
if stree[cls_id].nodeExist(event_id):
FOUND = True
break
if FOUND:
print 'root_id=%s'%root_id
else:
print 'it is a background event'
triggering cascade¶
In [79]:
#--- leaf vs. depth
import networkx as nx
#import treeDataStruct as tds
from itertools import count
import random as rnd
#dd={'ROOT':16005,
# 'exponent':4.0}
#dd={'ROOT':15724,
# 'exponent':4.0}
dd={'ROOT':9367,
'exponent':4.0}
#for ii in xrange(1):
ROOT = dd['ROOT']
exponent=dd['exponent']
for item in sl:
cls_id = item[2]
root_id = item[1]
if root_id != ROOT:
continue
tree.pairs = []
stree[cls_id].plot(ax)
tree.pairs = [root_id]+tree.pairs
for item in sl:
sid = item[2]
root_id = item[1]
if root_id != ROOT:
continue
print 'root_id=',root_id, 'n=%s'%(item[0])
tree.g_netx = None
tree.g_netx = nx.DiGraph(directed=True)
tree.g_netx.add_node(stree[ sid ].data,color='b')
stree[ sid ].sprint(swarm_lohi )
# read edgelist from grid.edgelist
nx.write_edgelist(tree.g_netx, path="grid.edgelist", delimiter=":")
H = nx.read_edgelist(path="grid.edgelist", delimiter=":") #,create_using=nx.DiGraph())
# nx.draw(H,with_labels=True,arrows = True)#pos=nx.spring_layout(H))
# nx.draw_networkx(tree.g_netx, arrows=True, with_labels=True,node_color='b') #, **kwds)
pos = nx.spring_layout(H)
nodes=H.nodes
colors = ['red' for n in nodes]
sizes=[exponent**swarm_lohi.loc[ int(n) ]['magnitude'] for n in nodes]
node_shape = ['o' for n in nodes]
for n,indx in zip(nodes,xrange(sys.maxint)):
if int(n) == root_id:
colors[indx] = 'yellow'
# node_shape[indx]='*'
fig=plt.figure(figsize=(8,8))
ax=fig.add_subplot(111)
ec = nx.draw_networkx_edges(H, pos, alpha=0.2,arrows = True, ax= ax,edge_color='b',width=3)
nc = nx.draw_networkx_nodes(H, pos,node_size=sizes, node_color=colors,cmap=plt.cm.seismic, ax=ax,
edgecolors='black',node_shape='o') #, nodelist=nodes, ,
# with_labels=False, node_size=100, cmap=plt.cm.jet)
plt.axis('on')
ax.tick_params(labelbottom=False,labelleft=False,length=0.0)
fig.show()
# DrawFrame(ax, (0.0,0.0),(0.0,0.0),.1)
# DrawFrame(ax, (0.16,0.16),(0.07,0.03),0.1)
DrawFrame(ax, (0.16,0.1),(0.09,0.03),0.1)
#plt.tight_layout()
#fig.savefig('%s/tree.png'%DIR_OUTPT_figs,dpi=300,bbox_inches='tight',transparent=True)
timeseries associated with the trig. cascade¶
In [80]:
#--- plot clusters
def Inside(t,(tmin,tmax)):
if tmin<= t<tmax:
return True
return False
def GetTrigList( nij_trig ):
d_trig = {}
for items in nij_trig.itertuples():
triggerID = items.Parent_id
d_trig.setdefault( triggerID, [] ).append( items.Event_id )
return d_trig
def PlotArrowsTime(df_complete,d_trig, (tmin,tmax), (xmin,xmax), (ymin,ymax) ):
tlist=[]
mlist=[]
xm = []
ym = []
for triggerID in d_trig: #--- mother events
t0 = df_complete['date'].iloc[ triggerID ] #--- time of mother event
m0 = df_complete['magnitude'].iloc[triggerID] #--- magnitude
if not ( Inside(t0,(tmin,tmax) ) and #--- within the given interval
Inside(swarm_lohi['x(km)'].iloc[triggerID],(xmin,xmax)) and
Inside(swarm_lohi['y(km)'].iloc[triggerID],(ymin,ymax)) ):
continue
# x0 = md.date2num( t0 ) #--- convert to time object
# y0 = m0 #--- magnitude
# tlist.append( t0 )
# mlist.append( m0 )
for daughter_id in d_trig[triggerID]: #--- children
t1 = df_complete['date'].iloc[ daughter_id ] #--- time
m1 = df_complete['magnitude'].iloc[ daughter_id ] #--- magnitude
if not ( Inside(t1,(tmin,tmax) ) and #--- within the given interval
Inside(swarm_lohi['x(km)'].iloc[daughter_id],(xmin,xmax)) and
Inside(swarm_lohi['y(km)'].iloc[daughter_id],(ymin,ymax)) ):
continue
# x1 = md.date2num( t1 ) #--- convert to time object
# y1 = m1 #--- magnitude
# xw = x1 - x0 #--- t_child-t_parent
# yw = y1 - y0 #--- m_child - m_parent
if daughter_id in tree.pairs:
tlist.append( t1 )
mlist.append( m1 )
ax.annotate("", (t1,m1), xytext=(t0, m0), zorder=1,
textcoords='data',clip_on=True,
arrowprops=dict(arrowstyle="-|>,head_width=%s,head_length=%s"%(head_width,head_length),
shrinkA=shrinkA,
color="b",linewidth="0.3"))
#--- plot mother
if triggerID == tree.pairs[0]: #in roott:
xm.append(t0)
ym.append(m0)
#--- plot circles
df=pd.DataFrame({'date':tlist,'mag':mlist})
ax.scatter(df['date'], df['mag'],
s=exponent**(df['mag']),
alpha=1,zorder=2,
facecolors='red',color='black')
##plt.savefig('timeSeries.png')
for x0,y0 in zip(xm,ym):
# print x0,y0
ax.scatter(x0,y0,
s=exponent**y0,
marker='*', zorder=3,
facecolors='yellow',color='black',
alpha=1.0)
#--------------------------------------
#----- key: event value: aftershock id
#--------------------------------------
d_trig = GetTrigList( nij_trig )
# #--------------------
# #----- subset (the same as trig. analysis)
# #--------------------
# swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
# 'magnitude',
# ( mc, sys.maxint ) )
# swarm_lohi.reset_index(inplace=True,drop=True)
# #--------------------
# #--- cartesian coords
# #--------------------
# swarm_lohi = GetCartesian( swarm_lohi )
# dd={'ROOT':16005,
# 'exponent':4.5,
# 'dt_before':0.06,'dt_after':1,
# 'head_width':.2,'head_length':.4,
# 'shrinkA':0}
# dd={'ROOT':15724,
# 'exponent':4.5,
# 'dt_before':0.06,'dt_after':1.5,
# 'head_width':.2,'head_length':.4,
# 'shrinkA':0}
dd={'ROOT':9367,
'exponent':4.5,
'dt_before':0.01,'dt_after':0.5,
'head_width':.2,'head_length':.4,
'shrinkA':0}
#for ii in xrange(1):
ROOT = dd['ROOT']
exponent=dd['exponent']
dt_before = dd['dt_before']
dt_after=dd['dt_after']
head_width = dd['head_width']
head_length=dd['head_length']
shrinkA=dd['shrinkA']
swarm_lohi = GetCartesian( swarm_lohi )
dff = swarm_lohi.sort_values('magnitude',ascending=False)
roott = [i[1] for i in sl]
#for item in dff.iloc[:30].itertuples(): #--- plot clusters associated with large events
for item in dff.itertuples(): #--- plot clusters associated with large events
if item.Index != ROOT:
continue
(xc,yc)=(item.longitude,item.latitude) #(-96.807159,35.985256) #(-98.739380,36.467762) #(-99.0,36.5)
tc = item.date #pd.to_datetime("2016-09-03 12:02:44.520")
dx = 1000
fig = plt.figure(figsize=(24,8),frameon=False) #, dpi=150)
ax = fig.add_subplot(111)
#--- plot all events
plt.ylim(mc,swarm_lohi['magnitude'].max())
plt.scatter(swarm_lohi['date'],swarm_lohi['magnitude'],
s=exponent**swarm_lohi['magnitude'],
alpha=0.1)
tt = tc #swarm_lohi['date'].loc[key_max]
xx={0:xc,1:yc}
PlotArrowsTime( swarm_lohi, d_trig,
( tt+datetime.timedelta(days=-dt_before),tt+datetime.timedelta(days=+dt_after)),
(xx[0]-dx,xx[0]+dx),
(xx[1]-dx,xx[1]+dx) )
ax.set_xlim(tt+datetime.timedelta(days=-dt_before),tt+datetime.timedelta(days=+dt_after))
ax.set_ylim(mc-0.1,6)
ax.tick_params(axis='y',labelsize=32,right=True)
ax.tick_params(axis='x',rotation=20,labelsize=28,top=True)
ax.xaxis.set_major_formatter(dates.DateFormatter('%d %b %Y \n %H:%M'))
#--- frame
DrawFrame(ax, (0.04,0.04),(0.22,0.04),0.1)
#fig.savefig('%s/timeSeries.png'%DIR_OUTPT_figs,bbox_inches='tight',dpi=150)
spatial map associated with the trig. cascade¶
In [81]:
#--- plot clusters
def Inside(t,(tmin,tmax)):
if tmin<= t<tmax:
return True
return False
def GetTrigList( nij_trig ):
d_trig = {}
for items in nij_trig.itertuples():
triggerID = items.Parent_id
d_trig.setdefault( triggerID, [] ).append( items.Event_id )
return d_trig
def PlotArrowsSpace(df_complete,d_trig, (tmin,tmax), (xmin,xmax), (ymin,ymax) ):
xlist=[]
ylist=[]
mlist=[]
xm=[]
ym=[]
m_mother=[]
for triggerID in d_trig: #--- mother events
t0 = df_complete['date'].iloc[ triggerID ] #--- time of mother event
( x0, y0 ) = ( swarm_lohi['longitude'].iloc[triggerID],
swarm_lohi['latitude'].iloc[triggerID] )
m0 = df_complete['magnitude'].iloc[triggerID] #--- magnitude
if not ( Inside(t0,(tmin,tmax) ) and #--- within the given interval
Inside(x0,(xmin,xmax)) and
Inside(y0,(ymin,ymax)) ):
continue
for daughter_id in d_trig[triggerID]: #--- children
t1 = df_complete['date'].iloc[ daughter_id ] #--- time
( x1, y1 ) = ( swarm_lohi['longitude'].iloc[ daughter_id ],
swarm_lohi['latitude'].iloc[ daughter_id ] )
m1 = df_complete['magnitude'].iloc[ daughter_id ] #--- magnitude
if not ( Inside(t1,(tmin,tmax) ) and #--- within the given interval
Inside(x1,(xmin,xmax)) and
Inside(y1,(ymin,ymax)) ):
continue
if daughter_id in tree.pairs:
xlist.append( x1 )
ylist.append( y1 )
mlist.append( m1 )
ax.annotate("", (x1,y1), xytext=(x0, y0), zorder=1,
textcoords='data', annotation_clip=False,
# arrowprops=dict(width=.2,headlength=.4,shrink=0.5,color='b'))
# backgroundcolor="w",
arrowprops=dict(arrowstyle="-|>,head_width=%s,head_length=%s"%(head_width,head_length),
shrinkA=shrinkA,
color="b",linewidth="0.3"))
#--- plot mother
if triggerID == tree.pairs[0]: # in roott:
xm.append(x0)
ym.append(y0)
m_mother.append(m0)
#--- plot circles
df=pd.DataFrame({'x(km)':xlist,'y(km)':ylist,'m':mlist})
ax.scatter(df['x(km)'], df['y(km)'],
s=exponent**(df['m']),
alpha=1, zorder=2,
facecolors='red',color='black')
for x0,y0,m0 in zip(xm,ym,m_mother):
plt.scatter(x0,y0,
s=exponent**m0,
marker='*',zorder=3,
facecolors='yellow',color='black',
alpha=1.0)
# #--------------------------------------
# #----- key: event value: aftershock id
# #--------------------------------------
# d_trig = GetTrigList( nij_trig )
# #--------------------
# #----- subset (the same as trig. analysis)
# #--------------------
# swarm_lohi = DataFrameSubSet( swarm, #--- complete catalog
# 'magnitude',
# ( mc, sys.maxint ) )
# swarm_lohi.reset_index(inplace=True,drop=True)
# #--------------------
# #--- cartesian coords
# #--------------------
# swarm_lohi = GetCartesian( swarm_lohi )
# dd={'ROOT':16005,
# 'exponent':4.0,
# 'dx':.0125 * 0.8,
# 'head_width':.15,'head_length':.3,
# 'shrinkA':0}
# dd={'ROOT':15724,
# 'exponent':4.0,
# 'dx':.0125 * 1,
# 'head_width':.15,'head_length':.3,
# 'shrinkA':0}
dd={'ROOT':9367,
'exponent':4.0,
'dx':.0125 * 1,
'head_width':.15,'head_length':.3,
'shrinkA':0}
#for ii in xrange(1):
ROOT = dd['ROOT']
exponent=dd['exponent']
dx = dd['dx']
head_width = dd['head_width']
head_length=dd['head_length']
shrinkA=dd['shrinkA']
#--- plot arrows
#for ii in xrange(1):
swarm_lohi = GetCartesian( swarm_lohi )
dff = swarm_lohi.sort_values('magnitude',ascending=False)
roott = [i[1] for i in sl]
#for item in dff.iloc[:30].itertuples(): #--- plot clusters associated with large events
for item in dff.itertuples(): #--- plot clusters associated with large events
if item.Index != ROOT:
continue
(xc,yc)=(item.longitude,item.latitude) #(-96.807159,35.985256) #(-98.739380,36.467762) #(-99.0,36.5)
tc = item.date #pd.to_datetime("2016-09-03 12:02:44.520")
dt = 1000 #--- setup
fig = plt.figure(figsize=(8,8)) #,frameon=False) #, dpi=300)
ax = fig.add_subplot(111)
#--- plot all events
ax.scatter(swarm_lohi['longitude'], swarm_lohi['latitude'], # 'x(km)'],swarm_lohi['y(km)'],
s=exponent**swarm_lohi['magnitude'],
alpha=0.1)
tt = tc #swarm_lohi['date'].loc[key_max]
xx={0:xc,1:yc}
PlotArrowsSpace( swarm_lohi, d_trig,
( tt+datetime.timedelta(days=-dt),tt+datetime.timedelta(days=+dt)),
(xx[0]-dx,xx[0]+dx),
(xx[1]-dx,xx[1]+dx) )
ax.set_xlim(xx[0]-dx,xx[0]+dx)
ax.set_ylim(xx[1]-dx,xx[1]+dx)
#--- set tick labels
xlist=np.linspace(xx[0]-dx,xx[0]+dx,4,endpoint=True)
ax.set_xticks(xlist)
xlist=[r'$%4.2f^{\circ}$'%i for i in xlist]
ax.set_xticklabels(xlist)
#
ylist=np.linspace(xx[1]-dx,xx[1]+dx,4,endpoint=True)
ax.set_yticks(ylist)
ylist=[r'$%4.2f^{\circ}$'%i for i in ylist]
ax.set_yticklabels(ylist)
ax.tick_params(labelsize=24)
ax.tick_params(axis='y',labelright=True,labelleft=False,right=True,left=True)
ax.tick_params(axis='x',top=True)#,left=True)
#--- frame
# DrawFrame(ax, (0.16,0.16),(0.07,0.03),0.1)
DrawFrame(ax, (0.1,0.16),(0.09,0.03),0.1)
#fig.savefig('%s/spatialMap.png'%DIR_OUTPT_figs,bbox_inches='tight',dpi=150)
In [82]:
#--- r vs t (including indirect aftershocks)
dd = {}
for cls_id in stree:
index_m = np.array(m).searchsorted(stree[cls_id].coord[2]) - 1 #--- group clusters based on magnitude range
dd.setdefault(index_m,[]).append(cls_id)
colors = ['black','red','green','blue','cyan','brown']
for group_m, col in zip(dd, colors ): #--- loop over magnitude ranges
#--- plot
fig = plt.figure(figsize=(4,4))
ax = fig.add_subplot(111)
ax.set_xlim(1,1e5)
ax.set_ylim(1e-3,1e2)
ax.set_xscale('log')
ax.set_yscale('log')
cls_ids = dd[ group_m ] #--- cluster ids for a given range
delta_t = []
delta_r = []
for cls_id in cls_ids: #--- loop over each cluster
tree.rtm_dic = {} #--- initialize distance-time list
stree[cls_id].get_coords() #--- update
t_root = stree[cls_id].coord[ 0 ] #--- relative to the root
r_root = np.array( stree[cls_id].coord[ 1 ] )
for key in tree.rtm_dic:
t_afs = tree.rtm_dic[key][0]
r_afs = np.array(tree.rtm_dic[key][1])
delta_t += [ ( t_afs - t_root ).total_seconds() ]
delta_r += [ ( ( ( r_afs - r_root ) ** 2 ).sum() ) ** 0.5 ]
ax.scatter(delta_t,delta_r,color=col,alpha=0.04)
plt.show()
